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Army  land  managers  and  environmental  plan¬ 
ners  must  estimate  runoff  and  sediment  yield 
from  small,  unpaged  watersheds  on  Army 
training  lands  to  assess  the  condition  of  the 
lands  and  to  evaluate  alternative  erosion  con¬ 
trol  plans.  The  U.S.  Army  Construction  En¬ 
gineering  Research  Laboratory  (USACERL) 
developed  the  Army  multiple  watershed  storm 
water  and  sediment  runoff  (ARMSED) 
simulation  model,  which  is  based  on  the 
MULTSED  model  and  has  been  adapted  for 
Army  use. 

ARMSED  is  a  single  event,  distributed,  deter¬ 
ministic  simulation  model  that  operates  on  MS- 
DOS  compatible  microcomputers  with  512  K 
RAM.  A  10-megabyte  hard  disk  is  recom¬ 
mended. 

This  report  provides  guidance  and  documenta¬ 
tion  for  ARMSED  users.  Volume  I  contains 
guidance  for  selecting  and  estimating  the  para¬ 
meters  and  values  needed  as  input  for  the 
model.  Volume  II  describes  the  various  sub¬ 
routines  of  ARMSED  and  the  modifications  that 
have  been  made  to  arrive  at  the  current  ver¬ 
sion.  As  the  user  base  for  ARMSED  expands, 
the  model  will  be  updated  and  modified  to 
incorporate  new  data  ideas. 
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ARMSED,  A  RUNOFF  AND  SEDIMENT  YIELD  MODEL  FOR 
ARMY  TRAINING  LAND  WATERSHED  MANAGEMENT 
VOLUME  It:  PROGRAM  DOCUMENTATION 


1  INTRODUCTION 
Background 

Army  land  managers  and  environmental  planners  must  estimate  runoff  and  sedi¬ 
ment  yield  from  small,  ungaged  watersheds  on  Army  training  lands.  These  estimates  are 
needed  to  help  assess  the  condition  of  the  lands  and  to  evaluate  alternative  erosion  con¬ 
trol  plans.  Because  estimating  runoff  and  sediment  yield  is  a  difficult  hydrologic  task, 
mathematical  computer  models  can  be  an  important  part  of  the  process.  The  U.S.  Army 
Construction  Engineering  Research  Laboratory  (USACERL)  developed  the  Army  multiple 
watershed  storm  water  and  sediment  runoff  (ARMSED)  simulation  model,  which  is  based 
on  the  MULTSED  model  developed  at  Colorado  State  University.  USACERL  conducted 
st  eies  of  MULTSED  to  test  the  formulation  and  sensitivity  of  the  model.1  ARMSED  is 
Army  tailored  adaptation  of  MULTSED. 

ARMSED  is  a  single  event,  distributed,  deterministic  simulation  model  that  oper¬ 
ates  on  MS-DOS  compatible  microcomputers  with  512K  RAM.  A  10-megabyte  hard  disk 
is  recommended. 

Objective 

This  report  provides  documentation  and  guidance  to  ARMSED  users.  Volume  II 
focuses  on  the  structure  of  the  current  model,  a  description  of  the  different  subroutines 
(processes),  modifications  made  to  the  model,  and  recognized  limitations  and  problems 
one  may  encounter  when  applying  the  model.  Volume  I  is  a  guide  for  estimating  the  input 
parameters  used  in  ARMSED. 

Approach 

This  volume  describes  the  current  state  of  the  ARMSED  model.  A  description  of 
the  model  structure  and  components  is  given  and  each  subroutine  is  described.  Modifica¬ 
tions  from  the  MULTSED  model  are  explained  and  limitations  are  described. 

Mode  of  Technology  Transfer 

The  ARMSED  program  is  available  on  a  5 i -in.  floppy  disk  and  can  be  obtained  by 
contacting  Mr.  Robert  E.  Riggins  at  USACERL-EN,  P.  O.  Box  4005,  Champaign,  IL 
61820-1305.  Telephone:  commercial  217/373-7234,  or  toll-free  800/USA-CERL  (outside 
Illinois),  800/252-7122  (within  Illinois).  ARMSED  will  be  fielded  under  the  Integrated 
Training  Area  Management  Program  as  part  of  the  Maintenance  and  Scheduling  Support 
System.  As  the  user  base  expands,  the  model  will  be  updated  and  modified  to  incorporate 
new  data  and  ideas. 


lH.  G.  Wenzel,  Jr.  and  C.  S.  Melching,  An  Evaluation  of  the  MULTSED  Simulation  Mod  . 
to  Predict  Sediment  Yield,  Technical  Report  N-87/27/ADA185615  (U.S.  Arn.  y 
Construction  Engineering  Research  Laboratory  [USACERL],  September  1987). 
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2  MODEL  STRUCTURE  AND  SUBMODELS 


Structure 


The  original  version  of  MULTSED  contained  three  submodels:  MSEDl,  MSED2,  and 
MSED3.  MSED1  processed  data  on  upland  watersheds  and  planes  and  MSED2  sorted/ 
ordered  the  output  from  MSED1  for  use  in  MSED3.  MSED3  processed  data  on  channel 
routing.  ARMSED  maintains  this  modular  structure,  but  M8ED2  and  MSED3  have  been 
merged  into  one  submodel  called  MSED3.  Current  listings  for  the  ARMSED  submodels 
are  in  Appendix  A. 


MSED1 

MSED1  is  the  primary  submodel  of  ARMSED.  The  basic  building  blocks  of  the  sub¬ 
model  are  referred  to  as  either  watershed  units  or  plane  units.  Watershed  units  consist 
of  a  single  channel  and  two  contributing  overland  flow  planes,  one  on  either  side  of  the 
channel.  Watershed  units  are  normally  located  at  the  upper  reaches  or  on  the  perimeter 
of  a  large  drainage  basin.  They  are  unique  because  the  planes  and  the  channel  all  have  a 
"no-flow  boundary"  condition  at  the  upstream  (upslope)  end.  Plane  units  that  are  part  of 
a  watershed  unit  can  be  modeled  separately.  Plane  units  not  associated  with  a  watershed 
unit  typically  contribute  to  a  channel  segment  that  does  not  have  a  no-flow  boundary 
condition.  The  channel  is  not  a  part  of  the  plane  and  must  be  analyzed  using  the  pro¬ 
cesses  in  MSED3.  Combinations  of  watershed  and  plane  units  are  used  to  represent  the 
geometry  and  other  characteristics  of  the  drainage  basin  under  study.  The  MSED1  sub¬ 
model  contains  a  main  program,  15  subroutines,  and  1  function  that  are  discussed  below. 

Program  MULTSED l 

This  is  the  main  program  for  the  MSED1  submodel.  In  this  program,  the  name  of 
the  input  file  is  requested,  some  indexing  information  is  read  from  the  data  file,  sub¬ 
rout  i  es  ANAWAT  and  OUT  are  called,  output  information  is  written  to  files  7  and  8  for 
later  use  in  MSED2,  and  some  conversions  are  made  to  the  data  for  input  and  output 
needs.  The  program  contains  labeled  common  block  SED.  The  output  information  con¬ 
tained  in  files  7  and  8  is  the  unit  id  =  ISEG  (variable  name  in  the  submodel),  sediment 
yields  =  YIELD  (for  each  sediment  size  fraction  =  NSED),  total  water  yield  =  QT,  and 
hydrograph  values  at  each  time  step  =  QDUM.  The  information  is  written  to  file  7  for 
plane  units,  and  to  file  8  for  watershed  units. 

Subroutine  ANAWAT 

This  is  the  primary  subroutine  of  the  MSED1  submodel.  In  this  subroutine,  INPUT, 
TEMP,  INTRCP,  CUTOFF,  FORWRD,  BACK,  TRANSP,  QLAT,  and  RESIST  are  called. 
The  subroutine  initializes  the  output  variables  and  converts  runoff  and  yield  values  for 
output.  The  rainfall  excess  used  in  routing  is  EXCES.  Labeled  common  blocks  are 
COVER,  DEPTH,  and  SED. 

Subroutine  INPUT 

This  subroutine  reads  the  remainder  of  the  input  file  that  was  initially  accessed  in 
Program  MULTSED1.  The  information  includes  all  of  the  geometric,  soils,  sediment, 
cover,  and  rainfall  data  for  the  response  units.  (The  structure  and  contents  of  the  input 
file  are  presented  in  Appendix  B.)  Data  on  hydraulic  conductivity  and  rainfall  rate  are 
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converted  to  inches  per  minute.  The  overland  flow  length  is  reduced  by  the  fraction  of 
depression  storage  on  the  plane.  The  representative  grain  sizes  to  be  used  in  the  sedi¬ 
ment  transport  computations  are  calculated  as  the  geometric  mean  between  the  succes¬ 
sive  sizes  of  the  grain  size  distribution  curve  that  was  read  into  the  model.  The  repre¬ 
sentative  fraction  is  calculated  as  the  difference  in  cumulative  fraction  between  succes¬ 
sive  sizes.  Labeled  common  blocks  are  COVER,  SOIL,  SED,  and  VAR1. 

Subroutine  TEMP 

This  subroutine  is  used  to  select  a  kinematic  viscosity  of  water  for  the  temperature 
that  was  read  into  the  submodel.  The  viscosity  is  then  used  to  adjust  the  hydraulic  con¬ 
ductivity  for  the  soil.  Labeled  common  block  is  SOIL. 

Subroutine  I NTRCP 

This  subroutine  is  used  to  subtract  cover  interception  from  the  rainfall  data  and 
readjust  the  rainfall  hyetograph  as  necessary.  If  interception  is  present,  the  duration  of 
the  rainfall  is  reduced  to  satisfy  the  interception  depth  and  not  the  intensity.  Therefore, 
portions,  or  all  of  the  beginning  hyetograph  increments  may  be  eliminated  to  satisfy 
intercept.  Labeled  common  block  is  COVER.  A  warning  is  printed  if  all  of  the  rainfall  is 
intercepted. 

Subroutine  CUTOFF 

This  subroutine  is  used  to  compute  the  rainfall  excess  array,  EXCES.  A  modified 
Green-Ampt  infiltration  equation  is  used,  which  is  called  as  function  DF.  The  subroutine 
computes  the  infiltration  on  a  time  step  or  incremental  basis.  For  time  periods  beyond 
the  last  rainfall  increment,  rainfall  excess  can  be  considered  equal  to  zero  or  equal  to 
some  infiltration  rate.  Currently,  the  rainfall  excess  is  set  to  zero  for  those  times.  If  an 
infiltration  rate  is  needed  to  better  describe  the  process  taking  place  (i.e.,  continued 
losses),  then  a  negative  rate  must  be  used  instead  of  zero.  For  most  practical  problems, 
this  rate  occurs  so  late  in  the  storm  that  using  a  value  other  than  zero  has  little  effect 
on  the  hydrograph  or  sediment  yield.  Labeled  common  blocks  are  SOIL  and  COVER. 

Function  DF 

This  function  is  accessed  from  subroutine  CUTOFF.  Given  the  depth  of  infiltrated 
water  from  the  previous  time  step  (hyetograph  increment)  and  soil  hydraulic  properties, 
the  function  returns  the  maximum  depth  of  water  that  could  be  infiltrated  in  the  current 
time  period.  When  passed  back  to  CUTOFF,  this  value  is  compared  with  the  rainfall 
depth  in  the  increment  and  the  lesser  of  the  two  is  the  amount  that  is  infiltrated  during 
the  time  period.  If  the  rainfall  depth  is  greater  than  the  infiltration  depth,  then  water 
(EXCES)  is  available  for  routing.  If  all  the  rainfall  excess  has  been  infiltrated,  a  warning 
message  is  printed. 

Subroutine  FORWRD 

This  subroutine  is  used  to  calculate  the  characteristics  downslope  (or  downstream, 
in  the  case  of  a  channel).  A  characteristic  can  be  thought  of  as  an  open-top  container 
that  moves  down  the  slope  collecting  rainfall  (excess)  until  it  reaches  the  downstream 
boundary  and  empties  the  water  it  has  collected.  There  are  an  infinite  number  of  char¬ 
acteristics  traveling  at  any  one  time,  so  calculating  values  for  each  one  is  not  practical. 
Because  the  rainfall  excess  is  represented  as  a  hyetograph,  each  characteristic  is 
assumed  to  begin  at  the  start  of  a  hyetograph  increment.  The  characteristic  is 


9 


propagated  in  a  piecewise  manner  across  the  plane  (channel)  for  each  time  increment.  If 
the  characteristic  goes  beyond  the  downstream  boundary  in  some  time  increment,  the 
correct  time  of  arrival  for  the  characteristic  is  calculated  from  the  last  computed  posi¬ 
tion  of  the  characteristic.  If  the  rainfall  excess  ends  before  a  characteristic  reaches  the 
downslope  boundary,  then  the  characteristic  is  extended  to  the  boundary  as  a  straight 
line  when  the  loss  rate  (refer  to  CUTOFF)  has  been  set  to  zero.  Otherwise  the  charac¬ 
teristic  is  propagated  with  a  negative  excess  (loss).  When  the  time  increment  reaches 
the  end  of  the  hyetograph,  subroutine  STOP  is  called  to  compute  the  arrival  times  of  the 
characteristics  that  originate  after  excess  ceases. 

Subroutine  STOP 

Similar  to  subroutine  FORWRD,  and  used  in  conjunction  with  FORWRD,  this  sub¬ 
routine  calculates  the  arrival  times  of  those  characteristics  that  begin  after  the  rainfall 
excess  ceases.  Those  characteristics  represent  the  final  bounding  characteristic  (zero 
value)  for  the  hydrograph  (i.e.,  no  more  flow). 

Subroutine  BACK 

This  subroutine  can  be  thought  of  as  the  reverse  of  subroutine  FORWRD.  In  this 
subroutine,  characteristics  begin  at  the  downslope  boundary  and  are  propagated  upslope 
to  intersect  the  no-flow  boundary.  This  propagation  is  necessary  because  the  character¬ 
istics  propagated  downslope  do  not  arrive  at  equal  time  intervals.  If  an  equally  timed 
incremental  hydrograph  is  desired,  characteristics  must  be  propagated  upslope  to  deter¬ 
mine  the  corresponding  starting  time.  The  starting  time  for  a  given  characteristic  is 
bounded  by  the  starting  times  previously  selected  in  subroutine  FORWRD.  The  subrou¬ 
tine  uses  a  second  order  Newton’s  method  to  quickly  estimate  the  appropriate  staring 
time.  (See  Eggert2  for  a  recent  presentation  of  the  technique.)  Subroutine  ITEGR  is 
called  to  determine  the  excess  between  the  new  beginning  times  for  the  upslope  propa¬ 
gated  (backshot)  characteristics.  Subroutine  INDX  is  called  to  determine  the  bounding 
times  on  the  previously  computed  FORWRD  propagated  characteristics.  Labeled  com¬ 
mon  block  is  DEPTH.  The  time  adjusted  discharge  values  are  contained  in  array  QOUT. 

Subroutine  ITEGR 

This  subroutine  is  used  to  find  the  rainfall  excess  (plane)  or  lateral  inflow  (channel) 
that  forms  the  upslope  (upstream)  boundary  condition  for  the  BACK  subroutine.  If  the 
backshot  characteristics  bound  a  change  in  the  excess  hyetograph  rate,  the  new  cor¬ 
responding  excess  is  time  weighted  (i.e.,  total  depth  is  used). 

Subroutine  INDX 

This  subroutine  is  used  to  find  the  bounding  characteristic  time  for  the  BACK  sub¬ 
routine.  Iteration  is  used  until  the  correct  time  is  found. 


2K.  G.  Eggert,  ’’Upstream  Calculation  of  Characteristics  for  Kinematic  Wave  Routing,” 
Journal  of  Hydraulic  Engineering,  Vol  113,  No.  6  (American  Society  of  Civil  Engineers 
[ASCE],  June  1987). 
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Subroutine  QLAT 

This  subroutine  is  only  used  when  a  channel  in  a  watershed  unit  is  routed.  Overland 
flow  from  the  contributing  planes  is  summed  for  each  time  increment.  The  contributing 
planes  had  previously  been  routed  and  adjusted  with  BACK.  The  lateral  inflow  then 
serves  the  same  purpose  as  the  rainfall  excess  in  subroutines  FORWRD  and  BACK.  The 
adjusted  lateral  inflows  are  in  array  QLT. 

Subroutine  RESIST 

This  subroutine  is  called  to  calculate  kinematic  wave  routing  parameters  using 
Manning's  equation.  The  original  MSED1  contained  both  Manning's  and  Chezy's  equa¬ 
tions.  Because  Chezy's  equation  is  so  rarely  used,  it  was  removed  from  the  current 
version. 

Subroutine  TRANSP 

This  subroutine  represents  a  major  consolidation  of  the  sediment  yield  estimation 
routines  for  the  plane  and  the  channel  units  that  were  in  the  original  MULTSED  version. 
When  transport  is  called  by  program  MULTSED1,  a  parameter  (ID)  is  passed  that  tells  the 
subroutine  to  determine  the  yield  using  the  hydraulics  and  other  characteristics  of  a 
plane  (ID=0  or  1)  or  of  a  channel  (ID=2).  The  subroutine  uses  the  Meyer-Peter,  Mueller 
bed  load  and  the  Einstein  suspended  load  equations  to  determine  the  sediment  transport 
capacity.  Sediment  supply  in  a  plane  comes  from  raindrop  splash,  loose  sediment,  and 
overland  flow  detachment.  For  a  channel,  sediment  supply  comes  from  the  contributing 
planes  and  detachment  by  channel  flow.  Sediment  yield  is  the  lesser  of  the  two;  supply 
or  transport  capacity.  Transport  and  yield  are  determined  for  each  sediment  size  frac¬ 
tion  specified  by  the  input  data.  The  Einstein  equation  has  not  been  modified  from  the 
original  version  as  suggested  by  Melching  and  Wenzel3  because  the  original  derivation 
for  the  model  was  evaluated  and  found  to  be  adequate.  Labeled  common  blocks  are 
DEPTH,  SED,  SOIL,  and  COVER. 

Subroutine  POWER 

This  subroutine  is  used  to  calculate  the  necessary  integrals  for  the  Einstein  sus¬ 
pended  load  equation.  The  integrals  can  be  expanded  and  integrated  as  an  infinite  series 
of  terms.  These  terms  are  iterated  until  the  error  between  two  successive  terms  reaches 
0.001. 


Subroutine  OUT 

This  subroutine  controls  printing  the  results  to  the  CRT  and  to  a  disk  file  called 
SCREENl.OUT.  Sediment  yields  are  in  array  YIELD  and  discharge  rates  are  in  array  Q. 


MSED2 

This  submodel  has  been  merged  with  MSED3  as  a  subroutine  because  the  only  time 
it  is  needed  is  when  MSED3  is  run.  It  is  described  in  the  following  section. 


3C.  S.  Melching  and  H.  G.  Wenzel,  "Calibrating  Procedure  and  Improvements  in 
MULTSED,"  Hydraulic  Engineering  Series  No.  38  (Department  of  Civil  Engineering, 
University  of  Illinois,  July  1985). 
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MSKD3 


This  submodel  is  used  to  determine  water  and  sediment  transport  through  a  channel 
system  (one  or  more  channels)  that  is  supplied  by  the  v/alershed  and  plane  units  and  any 
upstream  channels  or  small  reservoirs.  The  submodel  contains  several  subroutines  that 
are  similar,  if  not  identical,  to  those  used  In  MSED1.  There  are  a  total  of  18  subroutines 
in  this  submodel.  The  submodel  uses  a  finite  difference  numerical  solution  to  solve  the 
continuity  and  momentum  equations  for  water  flow  in  a  channel  that  has  a  nonzero  up¬ 
stream  inflow.  The  analytical  solution  used  in  MSEDl  is  not  appropriate  here  because  of 
this  nonzero  condition.  Sediment  in  the  channel  is  routed  uncoupled  from  the  water  (i.e., 
erosion  and  deposition  do  not  affect  the  water  flow)  using  the  finite  difference  form  of 
the  continuity  equation  for  sediment.  Each  channel  is  subdivided  Into  five  equal  longitu¬ 
dinal  sections  for  the  routing  parts  of  the  model.  The  following  sections  describe  the  dif¬ 
ferent  subroutines  in  more  detail. 

Program  MULTSED.1 

This  is  the  main  driving  program  for  the  MSED3  submodel.  This  program 
interactively  queries  the  user  for  the  name  of  the  data  file.  Subroutines  MSED2,  DATAl, 
TEMP,  IN1TL,  and  ROUTE  are  called  from  this  program.  Output  files  SCREEN3.0UT  ar.d 
MSEDll.OUT  are  opened.  SCREEN3.0UT  receives  the  CRT  display  output  and 
MSEDll.OUT  receives  detailed  information  about  each  channel  segment  that  was 
modeled.  The  CRT  and  SCREEN3.0UT  only  have  information  about  the  last  (final)  chan¬ 
nel  modeled.  Labeled  common  blocks  in  MSED3  have  names  of  the  form  VAR*,  where  * 
is  a  positive  integer.  The  labeled  common  blocks  in  program  MULTSED-3  are  VAR1, 
VAR2,  VAR6  and  VAR7. 

Subroutine  MSED2 

This  subroutine  is  used  to  reorder  and  reformat  the  output  created  by  MSED1.  It 
essentially  changes  a  1  x  n  vertical  file  into  a  4  x  (n/4)  horizontal  and  vertical  file.  It 
uses  file  7  to  create  files  3  and  4,  and  file  8  to  create  files  13  and  14.  File  7  was  origin¬ 
ally  the  output  file  for  the  plane  units  and  file  8  was  the  output  unit  for  the  watershed 
units. 

Subroutine  DATAl 

In  this  subroutine,  the  channel  indexes  and  characteristics  are  read.  The  internal 
file  number  is  designated  as  9.  The  subroutine  also  reads  information  from  files  3,  4,  13, 
and  14  that  were  created  by  subroutine  MSED2.  Files  3  and  14  are  left  open  to  be  read 
later  in  the  subroutine.  This  subroutine  also  recognizes  and  reads  data  that  describe 
small  reservoirs  (stock  ponds)  in  the  channel  system.  In  these  cases,  reservoir  data 
instead  of  channel  data  are  read.  Labeled  common  blocks  are  VAR  1  through  10,  and  15 
through  18.  Some  conversions  of  the  reservoir  data  are  performed,  and  the  representa¬ 
tive  grain  sizes  are  calculated. 

Subroutine  TEMP 

The  function  of  this  subroutine  is  identical  to  the  one  found  in  submodel  MSED1. 
However,  in  this  subroutine,  the  hydraulic  conductivity  of  the  stream  bottom  is  modi¬ 
fied.  The  labeled  common  blocks  are  VAR  1,  2,  and  5. 
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Subroutine  INITL 

This  subroutine  is  used  lu  initialize  variables  in  the  submodel.  Several  values  are 
set  to  zero,  some  values  have  units  changed,  and  the  fall  velocity  for  sediment  particles 
is  calculated  for  later  use  in  subroutine  TRANSP.  Labeled  common  blocks  are  VAR  1 
through  9,  and  17. 

Subroutine  ROUTE 

This  subroutine  is  the  primary  computational  center  for  the  submodel.  From  here, 
subroutines  DATA2,  RES,  UPLAT,  RESIS,  CHINL,  WROUT,  PERTG,  TRANSP,  SROUT, 
CEASE,  and  OUT  are  called.  As  information  is  passed  back  to  ROUTE  from  the  called 
subroutines,  it  is  either  passed  directly  into  the  next  subroutine  called  or  it  is  modified. 
Subroutine  ROUTE  does  summations  and  conversions  on  the  output  before  calling  sub¬ 
routine  OUT.  Labeled  common  blocks  are  VAR  1,  2,  7,  8,  9,  and  17. 

Subroutine  DATA2 

This  subroutine  is  used  to  read  information  from  files  3  and  14  that  were  created  by 
subroutine  MSED2.  These  files  may  be  read  at  each  time  step,  if  needed.  Based  on  the 
flow  volume  weighting,  sediment  yields  from  the  plane  and  watershed  units  are  converted 
to  volume  concentrations  at  each  time  step.  Labeled  common  blocks  are  VAR  1  and  15. 

Subroutine  RES 

This  subroutine  is  called  by  ROUTE  if  the  segment  is  identified  as  a  reservoir.  This 
subroutine  uses  simple  inflow,  outflow,  and  change  in  storage  to  determine  the  volume  of 
water  that  passes  through  the  reservoir.  It  also  uses  a  trapping  coefficient  (supplied  in 
the  input  data)  to  determine  the  amount  of  sediment  deposited  in  the  reservoir.  Sedi¬ 
ment  that  is  not  completely  trapped  is  allowed  to  go  downstream  with  the  water  into  the 
next  channel.  Using  the  model  to  analyze  the  effects  of  reservoirs  on  a  watershed  is 
cautioned  against  because  of  the  simplistic  nature  of  this  subroutine  and  some  continuity 
errors  that  were  noticed  when  reservoirs  were  used  in  other  model  applications.  Labeled 
common  blocks  are  VAR  1,  2,  7,  8,  9,  and  17.  This  subroutine  calls  subroutine  UPRES. 

Subroutine  UPRES 

Tiiis  subroutine  is  used  to  identify  those  upstream  channels  and  watershed  units 
that  directly  contribute  to  a  reservoir.  Once  identified,  all  water  and  sediment  inflows 
are  summed  to  arrive  at  a  toted  inflow  for  the  time  period.  Labeled  common  blocks  are 
VAR  1,  3,  7,  8,  15,  and  17. 

Subroutine  UPLAT 

This  subroutine  is  used  for  channels  and  is  similar  to  UPRES  except  it  also  identi¬ 
fies  lateral  inflows  to  the  channel.  Once  identified,  the  water  and  sediment  loadings  dur¬ 
ing  the  time  period  are  computed.  Labeled  common  blocks  are  VAR  1,  3,  7,  8,  9,  10,  and 
15. 

Subroutine  CHINL 

This  subroutine  is  used  to  simulate  infiltration  in  the  channel.  It  has  been  modified 
from  the  original  MULTSED  version  that  infiltrated  too  much  water  because  of  summa¬ 
tion  problems  in  the  algorithms.  Infiltration  is  calculated  using  a  modified  Green-Ampt 
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form iilnt ion  that  incorporates  the  previous  flow  depth  in  the  channel  as  part  of  the  head 
term.  If  infiltration  is  not  desired,  a  hydraulic  conductivity  of  0.0  should  be  used  in  the 
data  set.  Labeled  common  blocks  are  VAR  1,  2,  9,  and  10. 

Subroutine  RESIST 

This  subroutine  is  the  same  as  that  in  submodel  MSED1. 

Subroutine  WROUT 

This  subroutine  is  the  water  routing  component  of  the  submodel.  It  uses  an  uncon¬ 
ditionally  stable  finite  difference  routing  technique  to  move  water  from  section  to  sec¬ 
tion  in  the  channel.  The  technique  uses  linear  and  nonlinear  schemes  to  iterate  to  a  solu¬ 
tion  that  satisfies  a  convergence  criteria  related  to  the  flow  area.  Labeled  common 
blocks  are  VAR  1,  8,  9,  and  10. 

Subroutine  PERTG 

This  subroutine  is  used  to  determine  how  much  sediment  is  available  on  the  stream 
bed  for  erosion  during  each  time  step.  If  sediment  is  available,  then  size  fractions  are 
recalculated  to  reflect  the  composition  of  the  materials  that  are  left.  As  the  finer 
materials  are  washed  away,  the  coarser  particles  armor  the  bed  of  the  stream  and,  unless 
flow  increases,  erosion  and  sediment  transport  decrease.  Labeled  common  blocks  are 
VAR  1,  2,  7,  8,  and  9. 

Subroutine  TRANSP 

This  subroutine  is  similar  to  TRANSP  in  submodel  MSED1,  except  that  only  the 
transport  capacity  is  determined  and  not  the  supply  and  resultant  sediment  yield. 
Labeled  common  blocks  are  VAR  1  through  5,  7  through  10,  and  18.  This  subroutine  calls 
subroutine  POWER. 

Subroutine  POWER 

This  subroutine  performs  the  same  function  it  performs  in  submodel  MSED1. 
Subroutine  SROUT 

This  is  a  very  complex  subroutine  that  is  used  to  route  sediment  through  each  of 
the  five  equal-length  subdivisions.  The  channel  is  subdivided  to  aid  routing.  At  each 
channel  reach,  the  sediment  transport  rate  as  determined  by  TRANSP  is  compared  with 
the  available  sediment  from  the  upstream  reach  or  inflow,  the  lateral  inflow,  and  the 
material  left  on  the  streambed  from  the  previous  time  period.  Sediments  that  can  be 
transported  are  transported  to  the  next  reach  and  those  that  cannot  be  transported  are 
deposited  on  the  stream  bottom  to  await  the  next  discharge.  Eventually  the  finer  sizes 
are  removed  from  the  channel  and  the  larger  sizes  remain.  Labeled  common  blocks  are 
VAR  1,  2,  4,  5,  7  through  10,  and  18. 

Subroutine  CEASE 

This  subroutine  is  used  when  lateral  inflow  has  ended  so  that  sediment  is  not  im¬ 
mediately  stopped.  Labeled  common  blocks  are  VAR  1,  4,  and  7  through  10. 
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Subroutine  OUT 


This  subroutine  is  used  to  write  the  model  output  from  all  the  channels  to  file  11, 
MSEDll.OUT,  and  to  write  output  from  the  last  (final)  channel  to  the  CRT  and  file  6, 
SCREEN3.0UT.  Labeled  common  blocks  are  VAR  1,  2,  6,  and  17. 
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3  MODKL  MODIFICATIONS  AND  LIMITATIONS 


Modifications 

The  MULTSED  model  has  undergone  several  modifications  since  the  original  version 
was  developed.  A  significant  change  was  to  adapt  the  model  to  run  on  MS-DOS  compat¬ 
ible  computers.  The  model  is  no  longer  restricted  to  mainframe  or  minicomputer  appli¬ 
cations.  The  model  has  been  streamlined  further  by:  (1)  removing  Chezy's  C  from 
MSED1  and  MSED3,  (2)  assigning  constant  values  to  DPOW  and  TAOCK  in  MSED1,  (3) 
assigning  constant  values  to  AGB,  BEX,  and  DELTS  in  MSED3,  and  (4)  removing  the  con¬ 
nection  unit  inputs  from  MSED3.  Data  files  have  been  combined  so  there  is  only  one 
user-supplied  file  for  each  MSED  submodel  rather  than  the  two  separate  files  used  in  the 
original  sulmodels.  The  formats  and  contents  of  the  new  data  files  are  listed  in  Appen¬ 
dix  B.  The  TRANSP  subroutine  in  MSED1  was  rewritten  to  combine  overland  and  channel 
flow  sediment  yield  calculations  in  one  subroutine  instead  of  the  two,  almost  identical, 
subroutines  in  the  original  version.  On  the  plane,  infiltration  was  turned  off  once  rainfall 
excess  ended,  and  in  the  channel,  infiltration  can  be  turned  off  by  entering  a  zero  for 
hydraulic  conductivity.  Submodel  MSED2  was  rewritten  as  a  subroutine  for  submodel 
MSED3.  Data  files  were  combined  so  that  indices  proceed  characteristics  in  a  single 
file.  Various  logic  and  computation  errors  discovered  by  other  researchers  who  applied 
the  model  have  bee'i  incorporated  as  appropriate. 


Limitations 

The  model  still  has  some  limitations  and  shortcomings.  First,  the  user  is  restricted 
to  75  response  units  for  MSED1  with  200  time  steps.  The  restriction  for  MSED3  is  a  total 
of  35  channels  and  reservoirs  and  200  time  steps.  Second,  experience  indicates  that 
watershed  units  should  not  exceed  2  sq  mi  and  channel  lengths  should  not  exceed  1  mi. 
Similarly,  channel  lengths  in  MSED3  should  not  exceed  1  mi.  Although  the  model  will 
still  operate  if  a  single  channel  is  longer  than  1  mi  or  a  single  watershed  is  larger  than  2 
sq  mi,  the  resultant  hydrographs  and  sediment  yields  do  not  make  any  sense  in  terms  of 
timing.  Larger  areas  and  lengths  should  be  avoided.  Third,  the  time  increment  in 
minutes  and  the  final  time  (duration)  of  the  hydrograph  seem  to  affect  the  peak  flow 
rates.  There  is  no  logical  explanation  why  this  has  occurred  in  some  applications,  but 
experience  has  indicated  that  the  shorter  the  time  step  (down  to  1  minute,  minimum),  the 
more  reliable  the  answer.  Fourth,  the  time  step  of  the  rainfall  hyetograph  can  affect  the 
results.  If  a  long  hyetograph  time  step  is  used,  the  resultant  infiltration  and  subsequent 
excess  rainfall  rates  may  not  be  estimated  very  well.  Therefore,  short  time  steps  are 
recommended.  Fifth,  the  sediment  transport  equations  used  in  the  model  are  not  appro¬ 
priate  when  shear  stress  caused  by  steep  slopes  and/or  high  flow  rates  greatly  exceeds 
the  stress  needed  to  move  the  particles.  Sixth,  routing  in  the  channels  can  be  greatly 
affected  by  the  choice  of  channel  cross-sectional  parameters.  Trying  to  force  the  chan¬ 
nel  to  act  as  a  rectangle  with  no  relationship  between  flow  area  and  wetted  perimeter 
causes  errors  in  the  outflow  rates  so  that  they  exceed  the  rainfall  rates.  Seventh,  the 
model  will  stop  without  warning  if  a  long  final  time  is  selected  for  the  hydrograph  dura¬ 
tion. 
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APPENDIX  A: 

SOURCE  CODE  LISTINGS 


PROGRAM  MULTSED1 
C 

C  THIS  PROGRAM  WAS  ORIGINALLY  DEVELOPED  AT  COLORADO  STATE 
C  UNIVERSITY  IN  THE  LATE  1970’S.  IT  HAS  BEEN  MODIFIED  BY 
C  RESEARCHERS  AT  NEW  MEXICO  STATE  UNIVERSITY  AND  THE 
C  UNIVERSITY  OF  ILLINOIS.  TIM  WARD  AT  NEW  MEXICO  STATE 

C  UNIVERSITY  HAS  CONVERTED  IT  TO  RUN  ON  THE  IBM  PC  USING 

C  MS-FORTRAN. 

C  LAST  UPDATE:  JUNE  1987  BY  TJW 

C 

COMMON/ SED/  DCOEFF , DPOW , DOF , CHDOF , TAOCK , NSED , PS I ( 10 ) , DMB I ( 10 ) 
DIMENSION  QDUM(200) , ITYPE( 75 ) ,TITL£(20) , ISEG( 75 ) , YIELD( 10,3) , 
+  TYIELD( 3 ) , IPRINT( 75 ) 

CHARACTERS  FNAME1 

C 

C  -  INITIALIZE  VARIABLES. 

C  EXTERNAL  INDEX 

TUP  =  0.0 

TIN  =  0.0 

TOUT  =0.0 
TRVOL=0 . 0 
TVINTR  =0.0 
TAREA  =  0.0 
DO  101  1=1,200 
QDUM( I )=0 .0 
101  CONTINUE 

C  ASK  FOR  THE  NAMES  OF  THE  INPUT  FILES 

WRITE(*,9000) 

9000  F0RMAT(30(/),30X, 'WELCOME  TO  NULTSEDl' , /, 

*23X, 1  PC  VERSION  V87.06  BY  TIM  J.  WARD',10(/)) 

WRITE(*,9001) 

9001  FORMAT(10X, 'WHAT  IS  THE  NAME  OF  THE  MSED1  INPUT  FILE?’ 

READ(*,9010)  FNAME1 
9010  PORMAT(A) 

WRITE(*,9012)  FNAME1 

9012  FORMAT(//,10X,' OPENING  FILE  ',A20,'  AS  THE  INPUT  FILE.',/, 

1  10X, ’OPENING  SCREEN1 .OUT  AS  THE  SCREEN  OUTPUT  FILE.’,/, 

2  10X, 'OPENING  MSED7.DAT  AS  AN  OUTPUT  FILE  FOR  MSED3 . ' , / , 

3  10X, 'OPENING  MSED8.DAT  AS  AN  OUTPUT  FILE  FOR  MSED3. ' ,/) 

C 

0PEN(1,FILE=FNAME1, STATUS* 'OLD' ) 

OPEN( 6, FILE*' SCREEN1 .OUT' , STATUS- ' NEW ' ) 

OPEN ( 7 , FILE* ' MS ED 7 . DAT ' , STATUS* ' NL„  ) 

OPEN( 8 , FILE* ’ MSED8 . DAT ' , STATUS* ’ NEW ' ) 

START  READING  THE  DATA 
READ( 1 , 4000 ) ( TITLE (I), 1-1,20) 

4000  FORMAT(20A4) 

WRITE  (6, 4001)(TITLE(I), 1-1,20) 

4001  FORMAT ( 1H  ,'  TITLE:  ',20A4) 

WRITE  (*,4002 )( TITLE (I ) , 1*1 , 15 ) 

4002  FORMAT(30(/),35X,’MSED1  IS  NOW  RUNNING.',/, 

*10X, 'THE  WATERSHED  TITLE  IS:  ' ,/ , 15X, 15A4 , 15(/>) 


19 


C  -  READ  IN  THE  TIME  INCREMENT  AND  FINAL  TIME  FOR  THE 

C  -  HYDROCRAPH. 

READ(1,4003)DTIM,FTIM 
4003  FORMAT ( 2  F 10 . 0 ) 

C 

C  — -  READ  NUMBER  OF  PLANES  (NPL), 

C  —  NUMBER  OF  SMALL  WATERSHEDS  (NWS) 

READ( 1 , I000)NPL,NWS 
1000  FORMAT (3110) 

C  —  READ  IN  TYPE  ARRAY  TO  IDENTIFY  THE  TYPE  OF  UNITS. 

NUMP=NPL+NWS 
DO  104  1=1 ,NUMP 

READ(  1,1000)ISEC( I ) , ITYPE( I ) , IPRINT(  I ) 

104  CONTINUE 

C  -  CONVERT  TIMES  TO  SECONDS 

DTIM=DTIM't60 . 

FTIM=FTIM*60. 

C  -  CALCULATE  NUMBER  OF  INCREMENTS  IN  HYDROCRAPH. 

NUM=IFIX(FTIM/DTIM)+1 

C  WRITE  THE  RESULTS  TO  THE  SCRATCH  FILES  FOR  USE  IN  MSED3 
WRITE(8,6000)NWS,NUM 
WRITE( 7 ,6000) NPL, NUM 
6000  FORMAT ( 6 110) 

C 

IF  ( NUMP .EQ.l.AND.NPL.EQ.l)  THEN 
C  —IF  NO  CHANNELS  NEED  TO  BE  ROUTED  THEN  USE  ONLY 

C  —  THE  SMALL  WATERSHED  MODEL  (ANAWAT). 

CALL  ANAWAT ( 1 , DTIM, FTIM, QDUM , TVINTR , TAREA , TRVOL , YI ELD , TYIELD , 

+  QT , NUM ) 

C  —  CONVERT  TOTALS  FROM  CFS  TO  AC-FT  UNITS. 

TV INTR=TVINTR/ 43560. 

TAREA=TAREA/43560. 

TRVOL=TRVOL/43560 . 

ID  =  0 

CALL  OUT( I SEC( K1 ) ,YIELD( 1,3), QDUM, NUM, QT, TAREA, TVINTR, NSED, DTIM, 
+  DMBI,TYIELD(3), TRVOL, ID) 

STOP  1 
ENDIF 
C 

C  START  THE  LOOP  FOR  THE  PLANES  AND  THE  WATERSHEDS 
C 

DO  109  Kl=l ,NUMP 
WRITE(*,9002 )  K1 

9002  FORMAT(//,20X, ’MSED1  IS  COMPUTING  FOR  RESPONSE  UNIT  ’,13,/) 

DO  110  IK2=1,NUM 
QDUM( I K2 )=0 . 0 
110  CONTINUE 

T1=0. 

T2=0 . 

T3=0 . 

C 

K2  =  I TYPE ( ISEC(Kl)) 

ID  =  K2 

GOTO  ( 1 1 , 12 ) ,K2 
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C  -  THIS  IS  THE  ONE  PLANE  CASE. 

11  CALL  ANAWAT(ISEG(K1),DTIM,FTIM,QDUM,T1,T2,T3, 

+  YIELD, TYIELD,QT, NUM) 

WRITE(7 ,6000)  ISEG(K1 ) ,NSED 

WRITE ( 7 ,6005 ) (YIELD(K3 , 1 ) ,K3=1 ,NSED) 

WRITE (7 ,6005)QT 

WRITE(7,6005)(QDUM(K3),K3=2,NUM) 

IF(IPRINT(ISEG(K1)).LT.0)  GO  TO  109 
Tl=Tl/43560. 

T2=T2/43560. 

T3=T3/43560. 

CALL  OUT( ISEG(Kl) , YIELD, QDUM , NUM , QT , T2 , T1 , NSED , 

+  DTIM,DMBI ,TYIELD,T3 , ID) 

GO  TO  109 

C  -  THIS  IS  THE  SUBWATERSHED  CASE  (2  PLANES  &  1  CHANNEL) 

12  CALL  ANAWAT(ISEG(K1),DTIM,FTIM,QDUM,T1,T2,T3, 

+  YIELD,TYIELD,QT,NUM) 

WRITE(8,6000)  ISEG(K1 ) ,NSED 
WRITE(8 ,6005 )( YIELD(K3 ,3),K3=1,NSED) 

WRITE(8,6005)  QT 

WRITE(8, 6005 )(QDUM(K3),K3=2, NUM) 

6005  FORMAT (4G15 . 6 ) 

IF(IPRINT(K1).LT.0)  GO  TO  109 
Tl=Tl/43560. 

T2=T2/43560. 

T3=T3/43560 . 

CALL  0UT(ISEG(K1),YIELD( 1,3), QDUM, NUM, QT,T2,T1, NSED, DTIM,DMBI, 
+  TYIELD(3) ,T3 , ID) 

109  CONTINUE 

REWIND  8 
REWIND  7 
CLOSE(8 ) 

CLOSE (7 ) 

STOP  '  MSEDl  IS  FINISHED' 

END 

C 
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n  n 


c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


v  i c  *  Yr  ic  ic  ic i  c i  r  v 


f  If  If  1C  1C  If  If  If  If  If  If  If  If  If  If  1C  If  If 


iV  iV  iV  iV  iV  if  iV  iV  iV  Vc  i V  iV  i’f  i’c  iV 


SUBROUTINE  ANAWAT( I F I LE , DTIM , FTIM , QOUT , TVINTR , TAREA, TRVOL , 
+  YIELD, TYIELD,QT,NUM) 

THIS  IS  THE  PRIMARY  CALLINC  SUBROUTINE  THAT  CONTROLS 
ALL  OF  THE  OTHER  SUBROUTINES. 

-  PARAMETER  DEFINITIONS 

I  FILE  =  SUBUNIT  NUMBER. 

DTIM  =  TIME  INCREMENT  FOR  HYDROCRAPH. 

FTIM  =  FINAL  TIME  OF  HYDROCRAPH. 

QOUT  =  ARRAY  OF  THE  RESULTING  DISCHARGES. 

TVINTR  =  TOTAL  VOLUME  OF  INTERCEPTION  IN  FT-"'-3. 
TAREA  =  TOTAL  AREA  IN  FT** 2 . 

TRVOL  =  TOTAL  VOLUME  OF  RAINFALL  IN  FT**3. 


COMMON/COVER/  CRNCOV( 2  ) , CANCOV( 2 ) , VC( ? ) , VC( 2 ) , FIMP(2 ) , SLOOSE( 2 ) 
COMMON/DEPTH/  Y(200) 

COMMON/ SED/  DCOEFF , DPOW , DOF , CHDOF , TAOCK , NSED, PS  I ( 10 ) , DMBI ( 10 ) 
DIMENSION  PLNCTH( 3),SLOPE(3), RAI NOD( 200 ),RAINT(200), 

+  EXCES ( 200 ), EXCEST( 200 ), RAIN ( 200 ),OUTIME( 200 ),RT( 200), 

+  QLT( 200 ) ,QLTIME( 200) ,QL( 200,2 ) ,QOUT( 200 ) , TYIELD( 3 ) , 

+  YIELD( 10,3) 

C 

C  —  INITIALIZE  VARIABLES. 

CHECK=0 . 0 
TRVOL=0 . 0 
TAREA  =0.0 
TVINTR  =  0.0 
ICHECK  =  0 
ERR  =  0.0 
TYIELD( 1 )=0.0 
TYI ELD( 2 ) =0 . 0 
TYIELD( 3 )=0 .0 
DO  101  1=1,200 
QL( I , 1 )=0 .0 
QL( 1 , 2  )=0 . 0 
QOUT(I)  =0.0 
QLT(I)  =  0.0 
QLTIME(I)  =  0.0 
Y ( I )=0.0 
101  CONTINUE 

DO  313  1=1,10 
Y I ELD( I , 1 ) =0 . 0 
Y I ELD( 1 ,2  )=0.0 
YIELD(I,3)=0.0 
313  CONTINUE 

C  —  INPUT  ALL  DATA  FOR  ANAWAT  FROM  THE  FIRST  FILE. 

CALL  I NPUT( I  PLANE , RA I NOD , RA I NT , NRA IN , PLNGTH , S LOPE , 

+  T , ADW , XN, A1 , B 1  ) 

NR=NRAIN 

-  CALCULATE  VISCOSITY  AND  CORRECT  HYDRAULIC  CONDUCTIVITY 

—  FOR  TEMPERATURE. 
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CALL  TEMP (T, VI SCO,  IPLANE ) 

C  -  CALCULATE  TOTAL  INCHES  OF  RAINFALL. 

TRINCH=RAINOD( 1 )*RAINT( 1 ) 

DO  100  1=2 ,NRAIN 

TRINCH=RAINOD( I )*(RAINT( I )-RAINT(I-l ) )+TRINCH 
100  CONTINUE 

-  THE  FOLLOWING  LOOP  CALCULATES  THE  EXCESS  OF  EACH 

-  PLANE  AND  ROUTES  THE  EXCESS  ON  EACH  PLANE. 

DO  102  1=1, IPLANE 
QT  =  0.0 

C  -  CALCULATE  AREA  OF  PLANE  (AREA). 

AREA=PLNCTH( I )*PLNGTH( 3 ) 

TAREA=TAREA+AREA 

C  -  CALCULATE  RAINFALL  VOLUME. 

TRVOL=TRVOL+AREA*TRINCH/ 12 . 

C  -  INITIALIZE  RAIN  ARRAY  EQUAL  TO  RAINOD  ARRAY. 

NRAIN=NR 

DO  103  J=1 , NRAIN 
RAIN(J)=RAINOD(J) 

RT( J)=RAINT( J) 

103  CONTINUE 

DO  113  J=1 ,200 
EXCEST( J)=0.0 
EXCES( J)=0 .0 
113  CONTINUE 

C  -  CALCULATE  THE  INTERCEPTION 

CALL  INTRCP ( ERR , RAI N , RT , NRAIN , VINTR , I ) 
TVINTR=TVINTR+VINTR*AREA/ 12 . 

IF(ERR.EQ.O.O)  GO  TO  10 
ERR=0 . 0 
GO  TO  50 

C  -  CALCULATE  THE  INFILTRATION. 

10  CALL  CUTOFF(ERR,EXCES,EXCEST,NEX, RAIN, RT, NRAIN, 

♦  I ,FTIM) 

IF(ERR. EQ.O . 0 .OR. FIMP( I ) .GT.O .0 )  GO  TO  20 
ERR=0 . 0 

C  -  FOR  THE  CASE  WHEN  THE  TOTAL  RAINFALL  IS  COMPLETELY 

C  -  INTERCEPTED  AND/OR  INFILTRATED. 

50  IF ( IPLANE. EQ. 1 )  RETURN 
ICHECK=ICHECK+1 
IF(ICHECK.EQ.2)  RETURN 
DO  104  IT=1 ,200 
QL(IT,I )=0.0 

104  CONTINUE 

GO  TO  102 

C  -  CONVERT  UNITS  OF  EXCESS  RAINFALL  IN/MIN  TO  FT/SEC 

C  -  AND  UNITS  OF  TIME  FROM  MIN.  TO  SEC. 

20  NNEX=NEX- 1 

DO  105  J=2 ,NNEX 

EXCES( J)=(EXCES( J)*( 1 .-FIMP( I ) )+RAIN( J-l )*FIMP(l))/720. 
EXCEST( J )=EXCEST( J )*60 . 

105  CONTINUE 

EXCES(NEX)=EXCES(NEX)*( 1 .-FIMP( I ) ) / 7  20 . 

EXCEST ( NEX ) =EXCEST( NEX ) *6  0 . 
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0  ---  DI-FINK  B.  B  IS  USKD  IN  THE  EQUATION  Q=A*Y**B. 

C  —  B=3 . 0  FOR  THE  CASE  OF  OVERLAND  FLOW. 

B=3 . 0 

C  —  DEFINE  A.  A  IS  USED  IN  THE  EQUATION  Q=A*Y**B. 

C  —  FOR  OVERLAND  FLOW  A  IS  A  FUNCTION  OF  THE  GROUND 

C  —  COVER,  SLOPE,  AND  VISCOSITY  OF  WATER. 

R=100 .+( ADW-I00 . )*(GRNCOV( I ) / 100 . )**2 . 

G=32 .174 

A=8 . *S  LO  P  E ( I ) *C / ( R*VI SCO ) 

C  -  CALL  THE  FORWARD  ROUTING  PROCEDURE. 

CALL  FORWRD( A , B , PLNCTH ( I ) , EXCES , EXCEST , NEX ,OUTIME ,TSTOP , FTIM 
+,  CHECK) 

I F ( I  PLANE . EQ . 2 )  CO  TO  60 

CALL  BACK( A , B , PLNCTH( I ) , DTIM ,TSTOP ,OUTIME , EXCEST, EXCES , 

+  QOUT, NEX, CHECK) 

ID  =  0 

CALL  TRANSPOD,  I,  SLOPE,  PLNCTH,  NR,  RAINT,RAINOD,R,  QOUT,  QL<  1,1), 
+  DTIM , FTIM , VISCO , A, B , YIELD ,TYI ELD) 

QT=0 . 0 

DO  201  JJ=2 , NUM 

QT=QT+(QOUT(JJ)+QOUT(JJ-l))/2. 

201  CONTINUE 

QT=QT*PLNGTH( 3 )*DTIM/43560 . 

DO  202  JJ=1 ,NSED 

YIELD( JJ ,1)=YIELD(JJ,1 ) /PLNGTH ( 3 ) 

202  CONTINUE 

RETURN 

60  CALL  BACK(A,B,PLNGTH(I), DTIM, TSTOP ,OUTIME , EXCEST, 

+  EXCES, QL( 1,1), NEX, CHECK) 

ID  =  1 

CALL  TRANS P( I D , I , SLOPE , PLNGTH , NR , RAINT , RAINOD , R , QIUT ,QL( 1 , 1 ) , 

+  DTIM, FTIM, VISCO, A, B, YIELD, TYI ELD) 

102  CONTINUE 

C  —  TO  ROUTE  THE  CHANNEL  CALL  QLAT  TO  SUM  THE  LATERAL 

C  -  INFLOWS. 

NQ=I FI X( FTIM/ DTIM) 

CALL  QLAT(QL,QLT,QLTIME,DTIM,NQ) 

IF(Al.NE.O.O)  CO  TO  215 

C  -  DEFINE  B.  B  IS  USED  IN  THE  EQUATION  Q=A*Y**B. 

C  —  FOR  A  CHANNEL  WITH  A  TRIANGULAR  CROSS  SECTION 

C  -  B=1 . 25 . 

B1=0 . 5 

-  DEFINE  A.  A  IS  USED  IN  THE  EQUATION  Q=A*Y**B. 

—  FOR  A  CHANNEL  A  IS  A  FUNCTION  OF  THE  CHANNEL 
C  — -  GEOMETRY  (Al),  CHANNEL  SLOPE  (SL0PE(3)),  AND  THE 

C  —  DARCEY-WEISBACH  RESISTANCE  FACTOR  (RESIST). 

Al  =  ( ( 2 . / ( ( 1 . /SLOPE( 1 ) )+( 1 . / SLOPE (2 ) ) ) )**0 . 5 )* 

+  ( ( 1 . + ( 1 . /SLOPE( 1 )**2 ) )**0 . 5 

+  +( 1 .+( 1 . /SLOPE( 2 )**2 ) )**0 . 5 ) 

215  CONTINUE 

CALL  RESIST(XN,Al,Bl,A,B,SLOPE(3)) 

C  -  ROUTE  THE  CHANNEL. 

CHECK=1 . 

CALL  FORWRD ( A , B , PLNGTH (3) ,QLT, QLTIME , NQ , OUTIME , TSTOP , FTIM 
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4*  CHECK) 

CALL  BACK( A , 8 , PLNGTH( 3 ) , DTIM , TSTOP , OUTIME ,QLTIME ,QLT , 

+  QOUT, NQ, CHECK) 

ID  =  2 

CALL  TRANSP ( I D , I , SLOPE , PLNGTH , NR , RAINT , RAINOD , R , QOUT ,QL( 1,1), 
+  DTIM, FTIM,VISC0,A1,B1, YIELD ,TYIELD) 

QT=0 .0 

DO  205  JJ=2,NUM 

QT  =QT+ ( QOUT ( J  J- 1 ) +QOUT ( J J ) ) / 1 . 

205  CONTINUE 

QT=QT"'DTIM/43560 . 

109  CONTINUE 
RETURN 
END 
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C  *  *  *  *  *  *  -,v  *  *  *  *  *  * 

C 


*** 


******************************* 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SUBROUTINE  I NPUT( I  PLANE , RAINOD , RAINT , NRAIN, 

PLNCTH, SLOPE, T, ADW,XN,A1 ,B1 ) 

—  THIS  SUBROUTINE  READS  WATERSHED  DATA  FOR  THE 

—  SUBROUTINE  ANAWAT.  DATA  IS  READ  FROM  THE  FILE 

—  WTTH  UNIT  NUMBER  =  1. 

—  PARAMETER  DEFINITIONS. 

I  PLANE  =  INDICATOR  OF  WHETHER  1-PLANE  OR  2-PLANES 

AND  1  CHANNEL  IS  IN  THIS  PARTICULAR  WATER¬ 
SHED. 

RAINOD=ARRAY  OF  RAINFALL  INTENSITIES  FOR  THE  STORM. 
RAINT  =  ARRAY  OF  FINAL  TIMES  CORRESPONDING  TO  THE 
RAINFALL  INTENSITIES. 

NRAIN  =  NUMBER  OF  RAINFALL  INCREMENTS  (INTENSITIES). 
PLNCTH  =  ARRAY  OF  PLANE  &  CHANNEL  LENGTHS. 

SLOPE  =  ARRAY  OF  PLANE  &  CHANNEL  SLOPES. 

T  =  TEMPERATURE  OF  WATER. 


COMMON/COVER/  CRNCOV( 2 ) , CANCOV( 2 ) , VC( 2 ) , VC( 2 ) , FIMP( 2 ) , SLOOSE(2 ) 
COMMON/SOIL/  WET  K( 2 ) , POROS( 2 ) , SAVE( 2  ) , SW( 2 ) , SI ( 2 ) , 

+  PLAS 1(3) ,COHM( 3 ) 

COMMON/SED/  DCOEFF , DPOW , DOF , CHDOF ,TAOCK ,NSED , PS  I ( 10 ) , DMBI ( 10 ) 
COMMON / VAR 1 /  TITL 
DIMENSION  TITL( 3 ) 

DIMENSION  RAI NOD( 200 ) , RAINT( 200 ) , PLNGTH( 3 ) , SLOPE( 3 ) , 

P(ll),D(il) ,DPRES( 3 ) , PIMP( 2 ) 

C 

C  —  READ  NUMBER  OF  PLANES. 

READ( 1, 1000 )( TITL (I), 1=1, 3) , I PLANE 
1000  FORMAT( / / 3A4 ,110) 

C  -  FOR  EACH  PLANE  READ  HYDRAULIC  CONDUCTIVITY, 

C  -  SOIL  POROSITY,  INITIAL  SOIL  SATURATION, 

C  -  FINAL  SOIL  SATURATION  AND  AVERAGE  SUCTION 

DO  100  1=1,2 

READ( 1,2000 )WET  K( I ) , POROS ( I ) , SI ( I ) , SW( I ) , SAVE( I ) , PLASI (I  )  , 

+  COHM(I) 

2000  FORMAT( 7F10.0) 

C  -  CHANGE  HYDRAULIC  CONDUCTIVITY  TO  UNITS  OF  IN/MIN. 

WET  K( I )=  WET  K( I )/60. 

100  CONTINUE 

C  -  READ  IN  CANOPY  AND  GROUND  COVER  DATA  FOR  EACH  PLANE. 

C  —  PIMP  IS  THE  PERCENT  OF  IMPERVIOUS  COVER  AND 

C  —  SLOOSE  IS  THE  AMOUNT  OF  LOOSE  SOIL  (POUNDS)  AVAILABLE  IN 

C  —  ADDITION  TO  THAT  WHICH  WILL  BE  SPLASH  OR  FLOW  DETACHED. 

READ( 1 , 3000 ) ( CANCOV( I ) , VC( I ) , CRNCOV( I ) , VG( I ) , 

+  PIMP(I), SLOOSE ( I ) , 1  =  1 , 2  ) 

FIMP( 1 )=PIMP( 1 ) / 1 00 . 

FIMP(2 )=PIMP(2 )/100. 

3000  FORMAT(6F10.0) 

C  -  READ  IN  SLOPE  AND  LENGTH  OF  EACH  PLANE. 

C  -  DEPRES  IS  THE  FRACTION  OF  THE  PLANE 

C  —  IN  DEPRESSION  STORAGE. 
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READ( 1 , 4000 ) ( SL0PE( I ) , PLNGTH( I ) , DPRES ! I ) , 1=1 , 2 ) 

READ( 1 ,4000 )  SLOPE(3),PLNGTH(3) 

4000  FORMAT(3F10.0) 

C  -  DECREASE  THE  LENGTH  OF  THE  PLANE  BY  THE  FRACTION 

C  -  OF  THE  DEPRESSION  STORAGE. 

PLNGTH( 1 )=PLNGTH( 1)*( 1 .-DPRES ( 1 ) ) 

PLNGTH( 2 )=PLNGTH( 2 )*( 1 . -DPRES ( 2 ) ) 

C 

C  -  READ  IN  CHANNEL  AND  RAIN/WATER  INDICES. 

C  -  IN  THE  ORIGINAL  VERSION,  CHEZY'S  C  WAS  READ  HERE. 

READ(1,5000)NRAIN,T,XN,A1,B1,ADW,COHM(3) 

5000  FORMAT(I10 ,6F10 .0) 

C  -  READ  IN  STORM  RAINFALL  DATA. 

DU  200  1=1 ,NRAIN 

READ( 1 ,6000)  RAINOD(I) ,RAINT(I) 

6000  FORMAT! 2F10.0) 

C  -  CONVERT  INTENSITY  FROM  IN/HR  TO  IN/MIN. 

RAINOD( I )=RAINOD( I ) / 60 . 

200  CONTINUE 

C  -  READ  IN  SEDIMENT  TRANSPORT/EROSION  CHARACTERISTICS. 

C  —  DPOW  AND  TAOCK  ARE  NOW  SET  AS  CONSTANTS 

READ( 1,1010) DCOEFF , DOF , CHDOF , PLASI ! 3 ) , NSED 
1010  FORMAT(4F10. 0,110) 

DPOW  =2.0 
TAOCK  =  0.047 

-  READ  IN  THE  SEDIMENT  SIZE  DISTRIBUTION 

-  AND  CALCULATE  THE  GEOMETRIC  MEAN  SIZE 

-  AND  FRACTION  FOR  EACH  INCREMENT. 

DO  210  1=1, NSED 

READ!1,6000)D(I),P!I) 

IF  d.GT.l)  THEN 
IB  =  I  -  1 

PSl! IB )=ABS! P(I )-P( IB)  ) 
DMBl!lB)=(D!l)*D!lB))**.5/304.8 
ENDIF 

210  CONTINUE 
NSED=NSED-1 
RETURN 
END 
C 
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c 

SUBROUTINE  TEMP(T,VISCO, IPLANE) 

C 

C  —  THIS  SUBROUTINE  CORRECTS  THE  VISCOSITY  AND 

C  —  HYDRAULIC  CONDUCTIVITY  FOR  TEMPERATURE  VARIATIONS 

C  —  FROM  THE  ASSUMED  TEMPERATURE  OF  68  DEGREES  (F). 

C  -  PARAMETER  DF.FINTIONS . 

C  T  TEMPERATURE  IN  DECREES  F. 

C  VISCO  =  KINEMATIC  VISCOSITY  (FT**2/SEC) 

C  IPLANF  =  NUMBER  OF  PLANES. 

C 

COMMON/SOIL/  WET  K( 2 ) , POROS( 2  )  , SAVE ( 2 ),SW(2),SI(2), 

+  PLAS 1(3), COHM ( 3 ) 

DIMENSION  TE( 10) ,V( 10) 

DATA  TE/ 32. ,40. ,50. ,60. ,68.  ,80. ,90. ,  100. ,120. , 140./, 

+V/ 1.93, 1.66, I. 4 1,1. 22, 1.09, 0.930, 0.826, 0.739, 0.609, 0.5 14/ 

C 

C  -  CALCULATE  NEW  VISCOSITY  BY  INTERPOLATION. 

DO  100  1=1,10 

IF(TE(I).LT.T)  CO  TO  100 
FAC1=(T-TE( I — 1 ) ) / ( TE ( I )-TE( 1-1 ) ) 

VISCO=V( I-1)+FAC1*(V(I)-V(I-1)) 

CO  TO  10 

100  CONTINUE 

C  —  ADJUST  THE  HYDRAULIC  CONDUCTIVITY. 

10  FAC2=VISCO/ 1 .09 
DO  101  J=l, IPLANE 

WET  K( J )=WET  K( J ) / FAC2 

101  CONTINUE 
VISCO=VISCO*.OOOOl 
RETURN 

END 

C 
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>************************************************************* 


SUBROUTINE  INTRCP( ERR , RAIN , RAINT, NRAIN , VINTR, I PL) 

—  THIS  SUBROUTINE  DETERMINES  THE  VOLUME  OF 

-  INTERCEPTED  RAINFALL.  INTERCEPTION  DEPENDS 

-  ON  THE  PERCENTAGE  OF  THE  GROUND  THAT  IS 

-  COVERED  BY  (CANCOV)  AND  GROUND  (GRNCOV),  AND 

—  THEIR  RESPECTIVE  WATER  HOLDING  CAPACITIES (VC, VG) . 

-  TOTAL  INTERCEPTED  VOLUME  =  VINTR. 

-  PARAMETER  DEFINITIONS. 

ERR  =  ERROR  INDEX. 

RAIN  =  ARRAY  OF  RAINFALL  INTENSITIES. 

RAINT  =  ARRAY  OF  FINAL  TIMES  FOR  EACH  RAINFALL 
INTENSITY. 

NRAIN  =  NUMBER  OF  RAINFALL  INTENSITIES. 

VINTR  =  TOTAL  AMOUNT  OF  RAINFALL  INTERCEPTED  IN  INCHES. 
I PL  =  INDICATOR  OF  WHICH  PLANE  THE  EXCESS  IS 
BEING  CALCULATED. 

COMMON/COVER/  GRNCOV( 2 ) , CANCOV( 2 ) , VG( 2 ) , VC( 2 ) , FIMP( 2 ) , SLOOSE( 2 ) 
DIMENSION  R  OLD(2O0) ,RAIN(200) ,RAINT(200) ,RTOLD(200) 

-  THE  R  OLD  IS  AN  ARRAY  WHICH  TEMPORARILY  STORES  THE 

—  STORM  INTENSITIES.  OPERATIONS  ARE  DONE  ON  THE  R  OLD 

-  ARRAY  TO  OBTAIN  AN  ARRAY  WHICH  EQUALS  THE  TOTAL 

-  RAINFALL  MINUS  THE  RAINFALL  THAT  IS  LOST  BY 

-  INTERCEPTION. 

-  INITIALIZE  R  OLD  ARRAY  EQUAL  TO  RAIN  ARRAY. 

DO  100  1=1, NRAIN 
R  OLD(l )=  RAIN(I) 

RTOLD( I )=RAINT( I ) 

100  CONTINUE 

—  CALCULATE  TOTAL  INTERCEPTED  VOLUME  OF  RAINFALL. 

VGC=GRNCOV ( I PL ) *VG ( I PL ) / 1 00 . 

VCC=CANCO V ( I P L ) *VC ( I PL ) / 1 00 . 

VINTR  =VGC+VCC 
VRAIN  =  0. 

-  THE  RAINFALL  LOST  TO  INTERCEPTION  IS  SUBTRACTED 

-  FROM  THE  TOTAL  RAINFALL. 

DO  101  1=1, NRAIN 
IF(I.EQ.I)  GO  TO  10 

VRAIN  =VRAIN+R  OLD(I )*(RAINT(I )-RAINT(l-l ) ) 

GO  TO  11 

10  VRAIN  =  VRAIN+R  OLD( 1 )*RAINT( 1 ) 

11  IF(VRAIN.GT. VINTR)  GO  TO  12 

RAIN(I)  =  0. 

101  CONTINUE 
GO  TO  13 

12  DRV  =VRAIN-VINTR 
DT=DRV/RAIN( I ) 

IF(I.EQ.l)  GO  TO  14 
RAINT( 1-1 )=RAINT( I )-DT 
RETURN 
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14  RAIN( 1 )=0 .0 

RAINT( 1 )=RAINT( 1 )-DT 
GO  TO  15 

C  —  PRINT  WARNING  IF  TOTAL  RAINFALL  IS  LESS  THAN 

C  —  TOTAL  INTERCEPTED  VOLUME  OF  RAINFALL. 

13  ERR=1 .0 

WRITE(* , 1000 )VRAIN 
WRITE(6 , 1000 )VRAIN 

1000  FORMAT( / / '  THE  ENTIRE  VOLUME  OF  RAINFALL  ',P10.3, 

+  '  INCHES,'/'  HAS  BEEN  ABSORBED  BY  INTERCEPTION.’/) 
CO  TO  16 

C  —  RESET  RAIN  ARRAY  EQUAL  TO  R  OLD  ARRAY. 

15  CONTINUE 
NRAIN=NRAIN+1 

DO  102  1=2 , NRAIN 
RAIN( f )=R  OLD(I-I) 

RAINT( I )=RTOLD( 1-1 ) 

102  CONTINUE 

16  RETURN 
END 

C 
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SUBROUTINE  CUTOFF ( ERR , EXCES , EXCEST , NEX , RAIN , RAINT , NRAIN , 

+  IPL,FTIM) 

-  THIS  SUBROUTINE  CALCULATES  THE  EXCESS  RAINFALL. 

-  INFILTRATED  RAINFALL  BASED  ON  THE  GREEN- AMPT  EQUATION. 

-  THIS  IS  NOT  A  CONTINUED  INFILTRATION  MODEL  (INFILTRATION 

-  DOES  NOT  CONTINUE  AFTER  THE  END  OF  THE  RAINFALL). 

PARAMETER  DEFINITIONS. 

ERR  =  ERROR  INDEX. 

EXCES  =  ARRAY  CONTIANINC  EXCESS  RAINFALL 
INTENSITY  VALUES. 

EXCEST  =  ARRAY  CONTI AN I NG  THE  FINAL  TIMES  FOR 
EACH  EXCESS  INTERVAL. 

NEX  =  NUMBER  OF  EXCESS  INTERVALS. 

RAIN  =  ARRAY  OF  RAINFALL  INTENSITIES. 

RAINT  =  ARRAY  OF  FINAL  TIMES  FOR  EACH  RAINFALL 
INTENSITY. 

NRAIN  =  NUMBER  OF  RAINFALL  INTERVALS. 

I PL  =  INDICATOR  OF  WHICH  PLANE  THE  EXCESS 

IS  BEING  CALCULATED. 

FTIM  =  FINAL  TIME  OF  HYDROGRAPH. 


COMMON/SOIL/  WET  K(2) ,POROS(2) ,SAVE(2) ,SW(2) ,SI(2) , 

+  PLASI(3),COHM(3) 

COMMON/COVER/  CRNCOV( 2 ) , CANCOV( 2 ) , VG( 2 ) , VC( 2 ) , FIMP( 2 ) , SLOOSE( 2 ) 
DIMENSION  RAIN( 200 ) , RAINT( 200), EXCES ( 200 ) , EXCEST( 200 ) 

-  INITIALIZE  VARIABLES. 

ITP=0 
FOl  =  0. 

F 02  =  0. 

T  =  0. 

FTIMM=FTIM/60. 

EXCES ( 1 )  =  0. 

EXCEST(l)  =  0. 

GAMMA=SAVE(lPL)*POROS(IPL)*(SW( IPL)-SI(IPL)) 

C  -  CALCULATE  NUMBER  OF  EXCESS  INCREMENTS. 

NEX=NRAIN+1 

C  -  THE  FOLLOWING  LOOP  ITERATES  EXCESS  INCREMENTS. 

DO  105  1=1, NRAIN 

C  —  CALCULATE  THE  RAINFALL  TIME  INTERVAL — DTM. 

IF(I.EQ.I)  GO  TO  11 
DTM=RAINT( I )-RAINT( 1-1 ) 

GO  TO  12 

11  DTM=RAINT( 1 ) 

12  T=T+DTM 

C  -  CALCULATE  THE  POTENTIAL  INFILTRATED  VOLUME. 

DELF  =  DF(F01 ,WET  K(IPL) , GAMMA, DTM) 

C  -  COMPUTE  THE  POTENTIAL  AVERAGE  INFILTRATION  RATE. 

FD  =  DELF /DTM 
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c 

c 

c 

c 

c 

c 


25 

105 


C 

c 

c 

c 

c 

c 

c 
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10000 

+ 

+ 


c 


R=RAIN( I ) 

-  COMPARE  THE  RAINFALL  INTENSITY  AND  THE  AVERAGE 

—  POTENTIAL  INFILTRATION  RATE.  IF  THE  RAINFALL 
-  INTENSITY  IS  GREATER  THAN  THE  INFILTRAION  RATE 

—  THEN  CALCULATE  EXCESS.  IF  THE  RAINFALL  INTENSITY  IS 

-  LESS  THAN  OR  EQUAL  TO  THE  INFILTRATION  RATE  THEN 

-  THE  EXCESS  IS  ZERO. 

IF(FD.GT.R)  THEN 
FD1S  =R 

EXCES( 1+1 )=  0.0 

ELSE 

ITP=1 

FD1S=FD 

EXCES( 1  +  1 )=R-FD1S 
ENDIF 

EXCEST(I+1)  =T 
FOl  =  FOl  +  FDlS-'-DTM 
F02=F01 
CONTINUE 

EXCEST(NEX+1)=1.E20 

NOTE:  IN  THE  ORIGINAL  VERSION  THE  EXCESS  RATE  AT  THE  NEXT  STEP 
WAS  SET  TO  A  NEGATIVE  KW.  THIS  MAY  NOT  BE  THE  BEST  WAY 
TO  CONSIDER  EXTENDED  INFILTRATION.  IT  IS  NOW  SET  TO 
A  ZERO  INFILTRATION  RATE.  RE-COMMENT  THE  FOLLOWING  LINES 
IF  YOU  WISH  TO  CHANGE  IT. 

EXCES(NEX+1)  =  0.0 
EXC£S(NEX+1)  =  -WET  K(IPL) 

NEX=NEX+1 

-  PRINT  WARNING  IF  THERE  IS  NO  EXCESS. 

IF(FIMP(IPL) .GT.O.O.OR.ITP.NE.O)  RETURN 
WRITE(6, 10000) 

WRITE(*, 10000) 

FORMAT ( / '  CUTOFF  FINDS  NO  RAINFALL  EXCESS.'/ 

'  NO  ROUTING  WILL  BE  ATTEMPTED.  CONTROL  RETURNED  TO'/ 

'  ANAWAT.'/) 

ERR=2 .0 

RETURN 

END 
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C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 
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FUNCTION  DF(F,WK,HEAD,DT) 

-  EXPLICIT  APPROXIMATION  FOR  GREEN-AMPT  INFILTRATION 

-  MODEL  DETERMINES  THE  POTENTIAL  INFILTRATION 

-  VOLUME  DURING  TIME  INCREMENT. 

-  PARAMETER  DEFINITIONS. 

F  =  ACCUMULATED  INFILTRATED  VOLUME  IN  INCHES. 

WK  =  HYDRAULIC  CONDUCTIVITY. 

HEAD  =  AV.  SUCTION*POROSITY*(FINAL-INITIAL  SATURATION) 

DT  =  TIME  INCREMENT  IN  MINUTES. 

A=1 . 

B=2.*F-WK*DT 
C=2.*WK*(HEAD+F)*DT 
DF=(-B+SQRT(B*B+4.*A*C ) ) / ( 2 .*A) 

RETURN 

END 
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c 

c 
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SUBROUTINE  FORWRD(A,B,D,QIN,TIMEIN,N,OUTIME,TSTOP,FTIM, CHECK) 

THIS  SUBROUTINE  IS  USED  TO  PROPAGATE  THE  CHARACTERISTIC  CURVES 
IN  A  DOWNSLOPE  (DOWNSTREAM)  DIRECTION  BEGINNING  AT  THE  STARTING 
TIME  OF  EACH  EXCESS  RAINFALL  OR  LATERAL  INFLOW  INCREMENT.  THE 
METHOD  USES  AN  ITERATIVE  SOLUTION  AND  PROGRESSES  FOR  EACH 
TIME  STEP  UNTIL  THE  DOWNSLOPE  BOUNDARY  IS  INTERSECTED.  THE 
INTERSECTION  IS  THE  TIME  OF  ARRIVAL  OF  THAT  CHARACTERI STIC . 

KEY  VARIABLES  ARE: 

QIN  =  ARRRAY  OF  EXCESS  RATES 
D  =  TOTAL  LENGTH  OF  PLANE  OR  CHANNEL 

X  =  POSITION  OF  CHARACTERISTIC  AT  THE  END  OF  A  TIME  STEP 

DEPTH  =  DEPTH  OF  FLOW  AT  A  POINT  ON  THE  CHARACTERISTIC 
OUTIME  =  ARRAY  OF  ARRIVAL  TIMES  FOR  THE  CHARACTERISTICS, 

USED  IN  SUBROUTINE  BACK 

DIMENSION  QIN(200),TIMEIN(200) , OUTIME ( 200 ) 

DO  10  K=  1 , 200 
OUT I ME ( K)~0 .0 
10  CONTINUE 
I F INAL=N- 1 
I  IFINL=r FINAL-1 
DO  200  K=1 , I IFINL 
KK=K 
XPRE=0 . 

X=0 . 

DEPTHS. 

DX=0 . 

DO  100  I  =K ,  I  FINAL 

DT=TIMEIN(l+l)-TIMEIN(l) 

IF(QIN(I+1).EQ.0.0)  GO  TO  50 
I F ( I . EQ . I  FINAL)  CO  TO  80 
DEPTH=DEPTH+QIN(I+1)*DT 
I F( DEPTH . LE . 0 . 0  )  DEPTH=0.0 
DUMX=DEPTH-QIN( 1+1 )*DT 
IF(DUMX.LE.O.O)  DUMX=0.0 
DX=(A/QIN( 1+1 ) )*(DEPTH**B-DUMX**B) 

XPRE=X 

X=X+DX 
GO  TO  60 
CONTINUE 

BECAUSE  OF  TRUNCATION  ERRORS  IN  THE  PC  VERSION,  THE  FOLLOWING 
DEPTH  TRAPS  WERE  PUT  IN  THE  MODEL  TO  CATCH  SMALL  NEGATIVE  VALUES. 

IF(DEPTH.LE.O.O)  DEPTH  =  0.0 
DX= A*B* ( DEPTH- * ( B- 1 . ) )*DT 
XPRE=X 
X=X+DX 

60  IF(X.LT.D)  GO  TO  100 
X=XPRE 
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DEPTH=DEPTH-QIN( I+l )*DT 
IF(QIN(I+1) .LE.O.O)  GO  TO  70 
I F ( DEPTH .LE.O.O)  DEPTH=0.0 

DELT=(((D-X)*QIN(I+1)/A+DEPTH**B)**(1./B)-DEPTH)/QIN(I+1) 

OUTIME(K)=DELT+TIMEIN(l) 

IF(OUTIME(K) .lt.ftim)  GO  TO  200 
TSTOP=FTIM 
RETURN 
70  CONTINUE 

IF(DEPTH. LE.O.O)  DEPTH=0.0 

OUTIME ( K ) = ( D-X ) / ( A*B* ( DEPTH** ( B- 1 . ) ) )+TIMEIN(I) 
IF(OUTIME(K). LT.FTIM)  GO  TO  200 
TSTOP=FTIM 
RETURN 
80  CONTINUE 

IF(DEPTH. LE.O.O)  DEPTH=0 .0 

xmax=x-(a*depth**b)/qin(i+i) 

C  THIS  IS  THE  POINT  WHERE  WE  NEED  TO  SEE  IF  WE  HAVE  RUN  OUT 

C  OF  TIME  STEPS  BEFORE  WE  HIT  THE  DOWNSLOPE  BOUNDARY.  IF  SO, 

C  GO  TO  SUBROUTINE  STOP. 

IF(XMAX.LT.D)  GO  TO  300 
I F ( DEPTH . LE . 0 . 0 )DEPTH=0 . 0 

DELT=( ( (D-X)*QIN( I  +  l ) /A+DEPTH**B)**( 1 . /B)-DEPTH)/QIN( I  +  l ) 
OUTIME(K)=TIMEIN( I )+DELT 
IF(OUTIME(K). LT.FTIM)  GO  TO  200 
TSTOP=FTIM 
RETURN 
100  CONTINUE 

IF( CHECK. EQ. 0.0. AND. QIN(N).LT. 0.0)  GO  TO  200 
IF(DEPTH. LE.O.O)  GO  TO  199 

OUTIME(K)=TIMEIN( IFINAL)+(D-X)/ ( A*B*( DEPTH** ( B- 1 . ) ) ) 

IF(OUTIME(K). LT.FTIM)  GO  TO  200 

TSTOP-FTIM 

RETURN 

199  OUTIME(K)=1.E30 
TSTOP=FTIM 
RETURN 

200  CONTINUE 

IF ( CHECK. EQ. 1.0)  GO  TO  310 
IF(QIN(N) .EQ.0.0)  GO  TO  310 
KK= IFINAL 
300  CONTINUE 

CALL  STOP(A,B,D,QIN,TIMEIN,N, OUTIME, TSTOP,FTIM,KK, CHECK) 

RETURN 

310  OUTIME( IFINAL)=1 . E30 
TSTOP=FTIM 
RETURN 
END 


***^^****************************************************** 


SUBROUTINE  STOP ( A , B , D , QI N , TIMEIN , N , OUTIME , TSTOP , FTIM , K , CHECK ) 
C 
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C  III  IS  SUBROUTINE  COMPLEMENTS  THK  FORWR1)  SUBROUTINE  IN  THAT 

C  IT  SOLVES  FOR  THE  CHARACTERISTICS  THAT  ORIGINATE  AFTER  THE 
C  EXCESS  HYETOCRAPH  OR  LATERAL  INFLOW  HAS  ENDED. 

C 

DIMENSION  QIN(200),TIMEIN(200),OUTIME(200) 

IFINAL=N-1 

IXNT=K-1 

IF(K.GT.l)  CO  TO  200 

DEPTH=0 . 0 

DO  100  1=2 , IFINAL 

DEPTH=DEPTH+QIN( I )*(TIMEIN( I )-TIMEIN( 1-1 ) ) 

100  CONTINUE 

OUTIME( 1 )=TIMEIN( IFINAL)-DEPTH/QIN(N) 

TSTOP=OUTIME( 1 ) 

IF(OUTIME(l).GT.FTIM)  TSTOP=FTIM 
RETURN 

200  CONTINUE 

DTT=(T1MEIN(K)-TIMEIN( IXNT))/2. 

TIME1=TIMEIN( jxnt)+dtt 
I  IF  I NL= I  FINAL- 1 

I F ( CHECK . EQ . 1 . 0 )  I IFINL=IFINAL 
300  CONTINUE 
DEPTH=0 . 

X=0. 

DO  400  I=IXNT,IIFINL 
II  =  I 

IF(I.NE.IXNT)  GO  TO  410 
DT=TIMEIN(K)-TIME1 
GO  TO  420 

410  DT=TIMEIN(I+1)-TIMEIN(I) 

420  IF(QIN(I+1).EQ.0.0)  GO  TO  430 
DEPTH=DEPTH+QIN(I+1)*DT 
IF( DEPTH . LE . 0 . 0 )  DEPTH  =  0.0 
DUMX=DEPTH-QIN( 1+1 )*DT 
IF(DUMX.LE.O.O)  DUMX=0.0 
DX=(A/QIN( 1+1 ) )*(DEPTH**B-DUMX**B) 

X=X+DX 

IF(X.GE.D)  GO  TO  435 
GO  TO  400 
430  CONTINUE 

I F( DEPTH . LE . 0 . 0 )  DEPTH  =  0.0 
DX=A*B*( DEPTH** ( B-l . ) )*DT 
X=X+DX 

IF(X.GE.D)  GO  TO  435 
400  CONTINUE 

IF( CHECK. EQ. 1.0)  GO  TO  435 
I F ( DEPTH . LE . 0 . 0 )  DEPTH  =  0.0 
XMAX=X-(A*DEPTH**B)/QIN(N) 

IF(XMAX.LT.D)  GO  TO  440 
CON=(XMAX-D)/D 
IF(CON. LT. 0 .05  )  GO  TO  500 
435  DTT=DTT/2. 

TIME1=TIME1+DTT 
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I F ( CHECK . EQ . 0 . 0 )  GO  TO  300 
DEPTH=DEPTH-QIN ( 1 1 + 1 ) *DT 
IF( DEPTH. LE. 0.0)  DEPTH  =  0.0 
IF(QIN(II+1) . EQ.0.0)  GO  TO  438 

DELT=(((D-X)*QIN(II+1)/A+DEPTH**B)**(1./B)-DEPTH)/QIN(II+1) 
GO  TO  439 

438  CONTINUE 

I F ( DEPTH . LE . 0 . 0 )  DEPTH  =  0.0 

OUTIME ( K ) =( D-X ) / ( A*B* ( DEPTH** ( B- I . ) ) )+TIMEIN ( I I ) 
OUTIME(K)=TIMEIN( II )+DELT 

439  CONTINUE 

IFCOUTIME(K).LT.FTIM)  GO  TO  300 

TSTOP=FTIM 

GO  TO  550 

440  DTT=DTT/2. 

TIME1=TIME1-DTT 
GO  TO  300 

500  CONTINUE 

I F( DEPTH . LE . 0 . 0 )  DEPTH  =  0.0 

DELT=*(  ( ( D-X )*QIN(N ) /A+DEPTH**B )**(  1 .  /  B)-DEPTH ) /QIN(N) 

OUTIME(K)=TIMEIN( II+l )+DELT 

IF(II .EQ. IXNT)0UTIME(K)=TIME1+DELT 

TSTOP=OUTIME(K) 

IF(TSTOP.GT.FTIM)TSTOP=FTIM 
550  DUMQ-QIN(K) 

N=N+I 

DO  600  I=K,N 
DUM1=TIMEIN(I) 

TIMEIN( I )=TIME1 

TIME1=DUM1 

DUM2=QIN(I) 

QIN(I)=DUMQ 
DUMQ=DUM2 
600  CONTINUE 
RETURN 
END 
C 
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SUBROUTINE  BACK  ( AL , BET , LEN , DTIM ,TSTOP ,TL , T ,Q, QOUT , NQ , CHECK ) 

THIS  SUBROUTINE  USES  THE  SUBDIVIDED  SOLUTION  DOMAIN 
AS  SUPPLIED  BY  SUBROUTINE  FORWRD  TO  CALCULATE  THE 
TIME  OF  ORIGIN  OF  CHARACTERISTICS  CORRESPONDING  TO 
AN  ARBITRARILY  SELECTED  TIME  ON  THE  DOWNSTREAM 
BOUNDARY.  IN  SHORT,  BACK  CALCULATES  CHARACTERISTICS 
IN  THE  UPSTREAM  DIRECTION.  THIS  ALLOWS  THE  DISCHARGE 
TO  BE  KNOWN  AT  CONVENIENT  TIME  INTERVALS  THEREBY 
FACILITATING  THE  FORMATION  OF  THE  LATERAL  INFLOW  TO 
THE  CHANNEL  AND  ALLOWINC  GREATER  EASE  IN  INTERFACING 
THIS  WATER  ROUTINC  SIMULATION  WITH  OTHER  WATERSHED 
PROCESS  MODELS. 

SUBROUTINE  BACK  PROCEEDS  AS  FOLLOWS: 

IT  FIRST  DETERMINES  THE  LOCATION  OF  THE  TIME  OF  INTEREST 
ON  THE  DOWNSTREAM  BOUNDARY  WITH  RESPECT  TO  THE 
CHARATERISTICS  SUPPLIED  BY  SUBROUTINE  FORWRD.  NEXT, 

THE  ROUTINE  FORMULATES  F(TO)  AND  ITS  FIRST  TWO 
DERIVATIVES,  FINALLY,  THIS  RESULT  IS  ITERATED 
UNTIL  A  SUITABLY  ACCURATE  ESTIMATE  OF  THE  TIME  OF 
ORICIN  IS  CALCULATE  USING  A  SECOND  ORDER  NEWTONS  METHOD. 


PARAMETER  DEFINITIONS 
AL  =  A  IN  Q=A*Y**B. 

BET  =  B  IN  Q=A*Y**B. 

LEN  =  SLOPE  LENGTH  (FT). 

DTIM  =  TIME  INCREMENT  FOR  HYDROCRAPH  (SEC). 

FTIM  =  ENDING  TIME  OF  HYDROGRAPH  (SEC). 

TL  =  ARRAY  OF  FINAL  TIMES  FOR  THE  CHARACTERISTIC 
LINES  FOUND  IN  FORWRD  SUBROUTINE  (SEC). 

T  =  TIME  ARRAY  FOR  INFLOW  (SEC). 

Q  =  INFLOW  ARRAY  (FT/SEC  OR  CFS/FT). 

QOUT  =  OUTFLOW  DISCHARGES  (CFS/FT  OR  CFS). 

NQ  =  NUMBER  OF  INCREMENTS  IN  INFLOW. 

TSTOP  =  TIME  WHEN  RUNOFF  STOPS. 


COMMON  /DEPTH/  Y (200) 

DIMENSION  TL(200),T(200),Q(200), CUMQ( 200 ) , QOUT ( 200 ) 
REAL  LEN 

ONE  TIME  INITIALIZATION  OF  LOOP  PARAMETERS  AND  TIMES. 
EPQ  =  l.E  -  5 

IF  ( CHECK. EQ. 1. )  EPQ  =  l.E  -  5  *  LEN 
NITER  =  10 
EP  =  0.0001 
TIME  =  0. 

B 1  =  BET  -  1. 

B2  =  BET  -  2. 

B3  =  BET  -  3. 

Y ( I )  =  0.0 

CALCULATE  CUMLATIVE  INFLOW  ARRAY. 

TOT  =  Q( I )  *  T( 1 ) 
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CUMQ(l)  =  TOT 
TEST=0.0 

DO  100  I  =  2,200 

TOT  =  TOT  +  Q(I )  *  (T( I )  -  T( I  -  1)) 

CUMQ(I)  =  TOT 
Y(  I )  =  0.0 

100  CONTINUE 

I  =  1 

110  IF  (TL(I).NE.TL(I  +  D)  GO  TO  120 
1  =  1  +  1 
GO  TO  110 

THIS  LOOP  ITERATES  FINAL  TIMES  FROM  0.0  TO  FTIM.  FOR 
EACH  FINAL  TIME  AN  INITIAL  TIME  (TEST)  IS  CALCULATED. 

120  CONTINUE 

DO  300  J  =  2,200 

TPRE=TEST 

TIME  =  TIME  +  DTIM 
IF  (TIME.GT.TSTOP)  RETURN 

THIS  SECTION  IS  USED  TO  CALCULATE  THE  DISCHARGE  FOR 
ALL  TIMES  LESS  THAN  TL(1).  FOR  THESE  TIMES  IT  IS  NOT 
NECESSARY  GO  CALCULATE  THE  UPSTREAM  TIME  OF  THE  CHARACTERIS 
SINCE  ALL  OF  THESE  CHARACTERISTICS  BEGIN  AT  T=0. 

IF  (TIME.LE.TL(l))  THEN 
DUMYT  =  0.0 

CALL  ITEGR  (DUMYT, TIME, DEPTHX,Q,T,NQ) 

QOUT(J)  =  AL  *  DEPTHX  *  *  BET 

Y(J)  =  DEPTHX 
GO  TO  300 
ENDIF 

130  IK  =  I 

FIND  THE  BOUNDS  FOR  THE  CHARACTERISTIC  LINE. 

DO  140  IJ  =  IK, 200 

IF  (TIME.LE.TL(IJ  +  1))  GO  TO  150 
1  =  1  +  1 

140  CONTINUE 

150  L2  =  I 

IF  THE  LATERAL  INFLOW  INTENSITY  IS  EQUAL  TO  ZERO, 

IT  IS  NECESSARY  TO  SKIP  UP  TO  THE  NEXT  CHARACTERISTIC. 

IF  (Q(L2  +  D.NE.O.)  GO  TO  160 
1  =  1  +  1 
GO  TO  150 

160  CALL  INDX(TIME,L,NQ,T) 

M  =  L  -  L2  -  1 

C  THE  FIRST  GUESS  OF  THE  TIME  OF  ORIGIN  TO  3E  CALCULATED 

C  (TEST)  IS  THE  AVERAGE  OF  THE  TIME  OF  ORIGIN  OF 

C  THE  LOWER  BOUNDING  CHARACTERISTIC,  AND  THE  SMALLER  OF 

C  THE  TIME  ORIGIN  OF  THE  UPPER  BOUNDING  CHARACTERISTIC 

C  AND  TIME. 

TUP  =  T(I  +  1) 

TDN=T( I ) 

I F ( TPRE . GT . TDN ) TDN=TPRE 

IF  (TUP. GT. TIME)  TUP  =  TIME 
TEST  =  (TDN  +  TUP ) / 2 . 

C  THIS  IS  AN  ANALYTICAL  SOLUTION  USING  A  THIRD  ORDER 
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C  NEWTON  APPROXIMATION  TO  SOLVE  FOR  THE  CHARACTERISTIC  LINE. 

DO  281  J1  =  I, NITER 

FF  =  -  LEN/(AL  *  BET) 

FTO  =  FF 
TLAST  =  TEST 
IF  (M.EQ.O)  THEN 
TEST  =  TIME 

+  -  ((LEN/(AL  *  (Q(L)  *  *  Bl)))  *  *  (l./BET)) 

CO  TO  290 
ENDIF 

IF  (M.EQ.I)  CO  TO  190 

Cl  =  Q(L2  +  1)  *  T(L2  +  1)  -  Q(L2  +  1)  *  TEST 
FTO  =  FTO  +  (l./(Q(L2  +  1)  *  BET))  *  (Cl  *  *  BET) 
FDTO  =  -  (Cl  *  *  Bl) 

FD2TO  =  Q(L2  +  1)  *  Bl  *  (Cl  *  *  B2) 

DO  180  J2  =  2  ,  M 
JJ  =  L2  +  J2 

A  =  -  Q(JJ)  *  T( JJ  -  l)  +  CUMQ( JJ  -  1)  -  CUMQ(L2) 

+  +  Q(L?  +  1)  *  T(L2) 

IF  (Q(JJ).EQ.O.)  THEN 
C1=-Q(L2+1)*TEST+A 
CDT=T( JJ )-T(JJ~l) 

FTO=FTO+CDT* (Cl **B  1 ) 

FDTO=FDTO+  B 1*( -Q( L2+ 1 ) ) *CDT* ( C 1 **B2 ) 
FD2TO=FD2TO+Bl*B2*Q(L2+l)*Q(L2+l)*CDT*(Cl**B3) 

ELSE 

Cl  =  Q(JJ)  *  T(JJ)  -  Q(L2  +  1)  *  TEST  +  A 
C2  =  Q( JJ )  *  T( JJ  -  1)  -  Q(L2  +  1)  *  TEST  +  A 
FTO  =  FTO  +  ( 1 . / ( BET  *  Q(JJ)))  *  ((Cl  *  *  BET) 

+  -  ( C 2  *  *  BET)) 

FDTO  =  FDTO  +  (  -  Q( L2  +  1)/Q(JJ))  *  ((Cl  *  *  Bl) 

+  -  (C2*  *  Bl)) 

FD2TO  =  FD2TO  +  (Q(L2  +  1)  *  Q(L2  +  1)/Q(JJ))  *  Bl 
+  *  ((Cl*  *  B2 )  -  (C2  *  *  B2 ) ) 

ENDIF 

180  CONTINUE 

C 

190  A  =  -  Q(L)  *  T(L  -  1)  +  CUMQ(L  -  1)  -  CUMQ(L2)  +  Q(L2  +  1) 

+  *  T(L2) 

IF  (Q(L).NE.O.)  CO  TO  200 
Cl  =  A  -  Q(L2  +  1)  *  TEST 
CDT  =  TIME  -  T(L  -  1) 

FTO  =  FTO  +  CDT  *  (Cl  *  *  Bl) 

FDTO  =  FDTO  +  CDT  *  Bl  *  (  -  Q(L2  +  1 ) )  *  (Cl  *  *  B2) 
FD2TO  =  FD2TO  +  CDT  *  Bl  *  B2  *  Q(L2  +  1)  *  Q(L2  +  1) 

+  *  (Cl*  *  B3 ) 

CO  TO  220 

200  Cl  =  Q(L)  *  TIME  -  Q(L2  +  1)  *  TEST  +  A 

C2  =  Q(L)  *  T(L  -  1)  -  Q(L2  +  1)  *  TEST  +  A 

IF  (C1.LT.O.O.OR.C2.LT.O.O)  GO  TO  260 

FTO  =  FTO  +  ( 1 . / (Q(L)  *  BET))  *  ((Cl  *  *  BET) 

♦  -  (C2  *  *  BET)) 

FDTO  =  FDTO  +  (  -  Q(L2  +  1)/(Q(L) >)  *  ((Cl  *  *  Bl) 

+  -  (C2  *  *  Bl)) 
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FD2T0  =  FD2T0  +  ((Q(L2  +  1)  *  Q(L2  +  1))/Q(U)  *  B1 
+  *((C1  *  *B2 )  -  (C2  *  *  B2) ) 


TEST  FOR  CONVERGENCE.  IF  TEST  IS  SUCCESSFUL,  THE 
ITERATIONS  ARE  COMPLETE.  THE  PROGRAM  THEN  PROCEEDS  TO 
THE  NEXT  TIME  INCREMENT.  IF  NOT,  THE  PROGRAM  USES 
THE  SOLUTION  TO  A  TRUNCATED  TAYLOR’S  SERIES  TO  CALCULATE 
A  NEW  TRIAL  VALUE  OF  TEST. 

220  IF  (ABS(FTO/FF).LE.EP)  GO  TO  290 

B  =  (2.  *  FDTO/FD2TO)  -  (2.  *  TEST) 

C  =  (2./FD2TO)  *  (FTO  -  FDTO  *  TEST 

+  +  0.5  *  TEST  *  TEST  *  FD2TO) 

D  =  B  *  B  -  4.  *  C 

IF  (D.LT.O.)  GO  TO  230 
DTEST  =  0.5  *  SQRT(D) 

GO  TO  240 


THE  NEW  ESTIMATE  OF  TEST  MUST  FALL  INSIDE  THE  REGION 
C  OF  THE  SOLUTION  DOMAIN  BOUNDED  BY  THE  CHARACTERISTICS 

C  SELECTED  ABOVE.  IF  IT  DOES  NOT,  THEN  AN  ESTIMATE  MUST 

C  BE  REPICKED  INSIDE  THAT  REGION.  THIS  PROCEDURE  MAY 

C  RESULT  IN  THE  REJECTION  OF  THE  ESTIMATE  PROVIDED  BY 

C  THE  NEWTON )S  METHOD.  AND  THE  SUBSTITUTION  OF  A  NEW 

C  NEW  AVERAGE  VALUE  AS  THE  NEW  GUESS  OF  TEST. 

230  CONTINUE 

TEST  =  TEST  -  FTO/ FDTO 

IF  ( TEST. LT. TUP. AND. TEST. GE.TDN)  GO  TO  280 
GO  TO  250 

240  DTEST  =  0.5  *  (SQRT(B  *  B  -  4.  *  C)) 

IF  (ABS(FTO/FF ) .LE.EP)  GO  TO  290 

TEST  =  0.5  *  (  -  B)  -  DTEST 

IF  ( TEST. LT. TUP. AND. TEST. GT.TDN)  GO  TO  241 

TEST  =  0.5  *  (  -  B)  +  DTEST 

IF  ( TEST. LT. TUP. AND. TEST. GE.TDN)  GO  TO  280 

TEST  =  TLAST 

GO  TO  230 

241  TESTB=. 5*(-B)+DTEST 

I F ( TESTB . GT . TUP . OR . TESTB . LT . TDN )  GO  TO  280 
TESTC=TLAST-FTO / FDTO 
TEST 1=ABS ( TESTC-TEST ) 

TEST2=ABS ( TESTC-TESTB ) 

IF(TEST1.GT.TEST2)TEST=TESTB 
IF(TEST.LT. TUP. AND. TEST. CT. TDN)  GO  TO  280 
TEST=TLAST 
GO  TO  230 


250 

IF  (TEST. GE. TUP) 

GO  TO  270 

260 

TEST  =  (TLAST  ♦ 

TDN ) / 2 . 

GO  TO  280 

270 

TEST  =  (TLAST  + 

TUP ) / 2 . 

280 

CONTINUE 

281 

CONTINUE 

C  THE  VALUES  OF  DEPTH  (OR  CROSS-SECTIONAL  AREA  FOR  CHANNELS) 
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C  AND  DISCHARGE  ARE  NOW  CALCULATED.  A  TEST  IS  MADE  TO  CHECK 

C  IF  THE  DISCHARGE  IS  NEGLIGIBLY  SMALL  OR  IF  THE  PRESELECTED 

C  DURATION  OF  THE  HYDROGRAPH  HAS  BEEN  EXCEEDED.  IF  SO  THE 

C  ROUTINE  RETURNS  TO  ANAWAT. 

290  CONTINUE 

CALL  ITEGR  (TEST, TIME, DEPTHX, Q,T,NQ) 

QOUT(J)  =  AL  *  DEPTHX  *  *  BET 

Y(J)  =  DEPTHX 

MD  =  NJQ 

IF  (CHECK. GT. 0. )  NP  =  NQ  +  1 

IF  (TIME.GE.T(NQ  -  1 ) .AND.QOUT( J) .LE.EPQ)  RETURN 
300  CONTINUE 
RETURN 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


100 

55 


101 


c 

c 

40 


102 

c 

c 

c 


SUBROUTINE  ITEGR(T1 ,T2 ,ANS, QIN, TIMEIN, N) 

-  SUBROUTINE  INTEGR  INTEGRATES  THE  EXCESS 

-  RAINFALL  (OR  INFLOW  IN  THE  CASE  OF  A  CHANNEL) 

-  HISTOGRAM  BETWEEN  AN  INITIAL  TIME  (Tl) 

-  AND  A  FINAL  TIME  (T2). 

—  PARAMETER  DEFINITIONS. 

Tl  =  INITIAL  TIME  OF  CHARACTERISTIC. 

T2  =  FINAL  TIME  OF  CHARACTERISTIC. 

ANS  =  AREA  BETWEEN  Tl  AND  T2  OF  IN  ARRAY. 

QIN  =  ARRAY  BEING  INTEGRATED. 

TIMEIN  =  ARRAY  OF  FINAL  TIMES  CORRESPONDING  TO  THE 
QIN  ARRAY. 

N  =  NUMBER  OF  ELEMENTS  IN  QIN  AND  TIMEIN  ARRAYS. 

DIMENSION  QIN(200),TIMEIN(200) 

-  INITIALIZE  DUMMY  VARIABLE  WHICH  IS  USED 

-  TO  STORE  INTERMEDIATE  ANSWERS. 

CUM  =0. 

-  FIND  INITIAL  TIME  INCREMENT. 

DO  100  1=2, N 
II=I 

IF(T1 . LT. TIMEIN(I) )  GO  TO  55 
CONTINUE 
IXNT=II 

-  FIND  FINAL  TIME  INCREMENT. 

DO  101  I=IXNT,N 

IF(T2 .GE.TIMEIN(I ) )  GO  TO  101 
IFIN=I 
GO  TO  40 
CONTINUE 
IFIN=N 

-  FIND  AREA  INBETWEEN  INITIAL  AND  FINAL 

-  TIME  INCREMENTS. 

CONTINUE 

DO  102  J=IXNT, IFIN 
CUMPRE=CUM 

CUM=CUM+(TIMEIN( J)-TIMEIN( J-l ) )*QIN( J) 

IF(J.EQ.IFIN)  CUM=CUMPRE+QIN( J)*(T2-TIMEIN( J-l ) ) 

CONTINUE 

-  CORRECT  INTERMEDIATE  ANSWER  BY  SUBTRACTING 

-  THE  AREAS  WHICH  SHOULD  NOT  BE  INCLUDED. 

IF(Tl.NE.O.O)  CUM=CUM-(T1-TIMEIN( IXNT-1 ) )*QIN(IXNT) 

-  SET  VALUE  OF  FINAL  ANSWER. 


ANS=CUM 

RETURN 

END 


C 

C  ^MrtHr**-AdHH^**^******VMn>**^^**^******iHr********************** 

C 


43 


SUBROUTINE  INDX(TIME, L,NQ,T) 

C 

C  —  THIS  SUBROUTINE  LOCATES  WHICH  INFLOW  TIME  INCREMENT  (L) 

C  —  CONTAINS  A  GIVEN  TIME  (TIME). 

C  —  PARAMETER  DEFINITIONS. 

C  TIME  =  TIME  OF  INTEREST. 

C  L  INDEX  OF  T  ARRAY  CONTAINING  THE  TIME. 

C  NQ  NUMBER  OF  INFLOW  TIME  INCREMENTS. 

C  T  ARRAY  OF  INFLOW  TIME  INCREMENTS  (SEC). 

C 

DIMENSION  T(200 ) 

C 

DO  100  1=1, NQ 
L=I 

IF(TIME.LE.T(I))  GO  TO  10 
100  CONTINUE 
L=NQ+1 
10  RETURN 
END 
C 
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SUBROUTINE  QLAT ( QL , QLT , QLT I ME , DT I M , KF ) 

-  THIS  SUBROUTINE  USED  ONLY  IN  THE  2  PLANES  AND 

-  1  CHANNEL  CASE.  IT  TOTALS  THE  LATERAL  INFLOW 

-  INTO  THE  CHANNEL. 

-  PARAMETER  DEFINITIONS. 

QL  =  DOUBLE  DIMENSION  ARRAY  CONTAINING  THE 
OUTFLOWS  FROM  EACH  CHANNEL. 

QLT  =  ARRAY  OF  THE  TOTAL  LATERAL  INFLOW. 

QLT I ME  =  ARRAY  OF  THE  TIMES  FOR  THE  TOTAL  LATERAL 
INFLOW  (QLT  ARRAY). 

DTIM  =  TIME  INCREMENT  USED  (SEC). 

KF  =  NUMBER  OF  ELEMENTS  IN  TIME  ARRAY. 

DIMENSION  QL(200 ,2) ,QLT(200) ,QLTIME( 200) 

-  CALUCULATE  THE  TOTAL  AVERAGE  INFLOW  OVER  A  GIVEN  TIME 

-  PERIOD. 

QLT(1)=0.0 
QLTIME( 1)=0.0 
KFF=KF+I 
DO  105  1=2, KFF 

QLT( I )=(QL(1 , 1 )+QL( 1 ,2 )+QL( 1-1 , 1 )+QL( 1-1 , 2)  )/2. 

QLTIME( I )=QLTIME(I-1 )+DTIM 
105  CONTINUE 
KFF=KFF+1 
KF=KFF 

QLT(KFF )=0 .0 

QLTIME(KFF)=QLTIME(KFF-1)+1.E30 

RETURN 

END 


45 


c 

SUBROUTINE  RESI ST(XN , Al , B1 , ALP , BET, SLP ) 

C 

C  —  THIS  SUBROUTINE  CALCULATES  THE  PARAMETERS  A  AND 

C  —  B  IN  THE  EQUATION  Q=A*AREA**B. 

C 

C  — -  THIS  IS  FOR  THE  MANNINGS  RESISTANCE. 

BET=(5.-2.*BI)/3. 

ALP=( ( SLP*2 . 2 1 ) / ( XN*XN*A1** ( 4 ./ 3 .)))**. 5 
RETURN 
END 
C 
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SUBROUTINE  TRANSP( ID , I , SLOPE , PLNGTH , NRAIN , RT , RAIN , R ,Q , QL , 

♦  DTIM,FTIM,VISC0,A1,B1, YIELD, TYIELD) 

SUBROUTINE  TRANSP  CALCULATES  THE  SEDIMENT  YIELD  BY  SIZE 
FOR  THE  PLANE  UNITS  AND  THE  CHANNEL  UNITS  (ID=0,1  OR  2). 

-  FOR  THE  PLANE  UNITS  ID=0  (SINGLE)  OR  1 (WATERSHED) 

THE  DECIMAL  PERCENT  OF  GROUND  AND  CANOPY  COVER  IS  USED  TO 
ESTIMATE  THE  OVERLAND  FLOW  ROUGHNESS.  RAINDROP  SPLASH  AND 
OVERLAND  FLOW  ARE  USED  TO  ESTIMATE  THE  AMONT  OF  MATERIAL 
AVAILABLE  FOR  TRANSPORT  AND  THE  AMOUNT  THAT  CAN  BE  TRANS¬ 
PORTED.  RILL  AND  GULLY  DEVELOPMENT  ARE  NOT  CONSIDERED. 

-  FOR  THE  CHANNEL  UNITS  ID=2 

THIS  SUBROUTINE  CALCULATES  THE  AMOUNT  OF  SEDIMENT  GENERATED 
IN  THE  CHANNEL  UNIT  OF  A  WATERSHED  DUE  TO  THE  INFLOW  FROM 
FROM  THE  PLANES.  THE  TOTAL  SEDIMENT  INFLOW  FROM  THE  PLANES 
IS  CALCULATED  AS  TSDLAT.  THE  INFLOW  FOR  EACH  SIZE  FRACTION 
IS  DETERMINED  AS  SEDLAT.  THE  TOTAL  TRANSPORT  CAPACITY  FOR 
EACH  SIZE  (TSEDQ)  IS  INITIALIZED. 

COMMON/DEPTH/  Y(20 0) 

COMMON/SED/  DCOEFF , DPOW , DOF , CHDOF , TAOCK ,NSED , PSI ( 10 ) , DMBI ( 10 ) 
COMMON/SOIL/  WET  K(2) ,POROS(2) ,SAVE(2) ,SW(2) ,SI(2) , 

+  PLASI(3),COHM(3) 

COMMON/COVER/  GRNCOV( 2 ) , CANCOV( 2 ) , VG( 2 ) , VC ( 2 ) , FIMP( 2 ) , SLOOSE ( 2 ) 
DIMENSION  Q(200) , YIELD( 10,3), SEDLAT ( 10 ) , SEDQ( 10 ) , TSEDQ ( 10 ) , 

+  SLOPE(3),PLNGTH(3),RAIN(200),RT(200),TYIELD(3),QL(200) 

DO  SOME  PRE-COMPUTATIONS  BEFORE  YOU  START  THE  SEDIMENT  LOOPS 
IF(ID.LE.l)  THEN 
C  THIS  IS  FOR  THE  PLANE  UNITS 
CG  =  GRNCOV(l)/100. 

CC  =  CANCOV( I ) / 100. 

REPOSE  =0.84 

C  INITIALIZE  VARIABLES. 

FI  =  FIMP(I) 

TSDLAT  =0.0 

C  CALCULATE  THE  AMOUNT  OF  SEDIMENT  DETACHED  FROM  RAINFALL. 

ARl  =  l.-CG-FI+CG*FI 

AR2  =  AR1*(1-CC) 

DAREA  =  AR2*PLNGTH(3)*PLNGTH(I) 

WTFAC  =  DAREA*(1.-P0R0S(I))*2. 65*62. 4/12. 

RTPRE  =0.0 
DO  55  J=l, NRAIN 

DQ  =  DCOEFF*(RAIN( J)*60 . )**DP0W 
TSDLAT  =  TSDLAT  +  DQ*(RT( J)-RTPRE)/60. 

RTPRE  =  RT( J) 

55  CONTINUE 

TSDLAT  =  TSDLAT* WTFAC  +  SLOOSE(I) 

ELSEIF  (ID.EQ.2)  THEN 
C  THIS  IS  FOR  THE  CHANNEL  CASE 
TSDLAT=0 .0 
DO  65  ISED=1,NSED 
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SEDLAT(ISED)=YIELD(ISED, 1 )+YI ELD( I SED , 2 ) 
TSDLAT=TSDLAT+SEDLAT( ISED) 

65  CONTINUE 
END  IF 

SSEDQ=0 .0 
DCOH=0 .0 
T=0 . 0 

DMBMAX  =  -1000.0 
DO  75  I J=1 ,NSED 
SEDQ( I J  )=0 . 0 
TSEDQ( I J  )=0 . 0 

IF(DMBKIJ)  .CE.  DMBMAX)  DMBMAX  =  DMBI(IJ) 

75  CONTINUE 

C 

C  LOOP  110  CALCULATES  THE  TRANSPORT  CAPACITY  FOR  ALL  THE  SIZE 

C  FRACTIONS  FOR  EACH  CALCULATED  DISCHARGE. 

DO  110  I J=2 ,200 

T=T+DTIM 

IF(T.GT.FTIM)  CO  TO  121 
IFOD.EQ.0.0R.ID.EQ.2)  QQ  =  Q(IJ) 

IF(ID.EQ.l)  QQ  =  QL(IJ) 

IF  (QQ.EQ.0.0)  CO  TO  110 
IF(ID.LE.l)  THEN 

C  CALCULATE  THE  DEPTH  OF  FLOW  AND  MEAN  VELOCITY 

DEPTHP  =  (QQ'--R*VISCO/(257.6*SLOPE(l)))**(l./3. ) 

VMEAN  =  QQ/DEPTHP 
C  EFFECTIVE  SHEAR  STRESS: 

TAO= . 24225*45 .*VISCO*QQ/ (DEPTHP**2 . ) 
IF(PLASKI).CT.O.O)  THEN 

TCOH=0 . 0022'-'-' PLASI  ( I  )**0 . 82 
ERATE=COHM( I )*(TAO/TCOH-l .  ) 

IF(ERATE.LE.O.O)  ERATE  =  0.0 
DCOH=DCOH+ERATE 
ENDIF 

SUMPS=0.0 

SK=COS(ATAN( SLOPE ( I ) ) )*SQRT( 1 .-(SLOPE ( I ) /REPOSE )**2 ) 
ELSEIF  (ID.EQ.2)  THEN 

C  THIS  SECTION  CALCULATES  THE  CHANNEL  CROSS- 

C  SECTIONAL  AREA  AND  TOP  WIDTH. 

C  THE  CROSS-SECTIONAL  AREA  IS  TRANSFERRED  AS  Y(IJ) 

CHAR=Y( I J) 

WPER=A1*CHAR**B1 
HYRAD=CHAR/WPER 
VMEAN =QQ/ CHAR 
TAO=62 .4*SL0PE( I )*HYRAD 
IF( PLASI ( I ) .GT.O .0 )  THEN 
TCOH=0 . 0022 ''-PLASI ( I  )** .  82 
ERATE=COHM( I )*(TAO/TCOH-l . )*WPER 
I F( ERATE.LT. 0.0)  ERATE=0.0 
DCOH=DCOH+ERATE 
ENDIF 

SUMPS  =  0. 

FMAX=  1.69+2.  -'-ALOG  10(2. '-'-HYRAD/ DMBMAX ) 

IF(FMAX.EQ.O.O)  FMAX=0. 0000001 
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c 

c 

c 


c 

c 


c 


c 


c 

c 

106 

c 


c 


c 

c 


FMAX=1.5/(FMAX)**2 
IF(FMAX.GT.O. 10)FMAX=0. 10 
IF(FMAX.LT. 0.010)  FMAX=0.010 
FMIN  =  FMAX/2. 

ENDIF 

THE  111  LOOP  CALCULATES  THE  BED  LOAD  AND  SUSPENDED 
LOAD  TRANSPORT  FOR  EACH  SIZE  FRACTION. 

DO  111  ISED=1 ,NSED 
DMB=DMBI(ISED) 

PS=PSI(ISED) 

CRITICAL  SHEAR  STRESS: 

IF(ID.LE.l)  THEN 
THIS  IS  FOR  THE  PLANE  CASE 

TAOC=102 .96*TAOCK*DMB*SK 
IF  (TAO.GT.TAOC)  GO  TO  106 
SEDQ(ISED)=0. 

SUSP=0. 

GO  TO  108 

ELSEIF  (ID.EQ.2)  THEN 
THIS  IS  FOR  THE  CHANNEL  CASE 

FDIVR=1.69+2.*ALOG10(2.*HYRAD/DMB) 

IF(FDIVR.EQ.O.O)  FDIVR=0 . 0000001 
F=  1 . 5  /  (FDIVR)*V*2 
I F ( F . GT . FMAX ) F=FMAX 
IF(F.LT.FMIN)  F=FMIN 
TAO= . 2425*F*VMEAN*VMEAN 
TAOC=102 .96*TAOCK*DMB 
CHECK  TO  SEE  IF  EXCESS  SHEAR  EXISTS 
IF  (TAO.GT.TAOC)  GO  TO  106 
SEDQ(ISED)=0. 

SUSP  =  0. 

GO  TO  108 

ENDIF 

BED  LOAD  CALCULATION  USING  THE  MEYER-PETER, MUELLER  EQUATION 
SEDQ(ISED)=( 12 .85/1 . 392)*(TAO-TAOC)**( 1 .5 ) 

SHEAR  VELOCITY: 

IF(ID.LE.l)  THEN 

SV=SQRT( 62 .4*DEPTHP*SLOPE( I ) / 1 . 938 ) 

ELSEIF  (ID.EQ.2)  THEN 

SV=SQRT( 62 . 4*HYRAD*SL0PE( 3 ) / 1 . 938 ) 

ENDIF 

FALL  VELOCITY  CALCULATION 

IF(DMB.GT. 0.0002)  THEN 

FVB=( SQRT( 36 . 064*DMB**3 . +36 . *VISCO**2 . )-6 . *VISCO)/DMB 
ELSF 

FVB=2 .951 7*DMB**2 . /VI SCO 
ENDIF 

CALCULATE  PARAMETERS  FOR  EINSTIEN'S  SUSPENDED 
SEDIMENT  LOAD  APPROXIMATION. 

ZR=FVB/(0.4*SV) 

BMV=2 . 5+VMEAN / SV 
IF(ID.LE.l)  THEN 
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AR=2  ."'DMBMAX/DEPTHP 
IF(AR.LT. .04)  AR= . 04 
IF  (ZR.GT.5.5.0R.AR.GT.0.5)  GO  TO  107 
ELSEIF  (ID.EQ.2)  THEN 
AR=2 . " DMBMAX/HYRAD 
IF  (AR.GT.0.9)  GO  TO  107 
ENDIF 

C  EVALUATE  J1  AND  J2. 

CALL  POWER  (ZR,AR,FJ,SJ,1.0E-2) 

P=AR  ZR-1 . ) / ( 1 1 . 6*( 1 . -AR )**ZR ) 

S  US  P = P* ( BMV " F  J + 2 . 5*S J ) 

SUPMAX= 1 . / AR 

IF(SUSP.GT.SUPMAX)  SUSP=SUPMAX 
IF  (SUSP. LE. 0.0)  SUSP=0.0 
GO  TO  108 

107  SUSP=0. 

108  CONTINUE 
C 

C  CALCULATE  THE  TOTAL  SEDIMENT  TRANSPORT  CAPACITY  FOR  THIS 

C  SIZE  FRACT I ON ( TS EDQ ) .  UPDATE  TRANSPORTING  CAPACITY  FOR 

C  FOR  THIS  SIZE  FRACTION  AT  THIS  TIME  INCREMENT  FOR  SU3PEN- 

C  DED  SEDIMENT.  UPDATE  THE  TOTAL  TRANSPORT  RATE  FOR  THIS 

C  TIME  INCREMENT. 

IF(ID.LE.l)  THEN 
WX  =  1.0  *  AR1 
PLX  =  PLNGTH( 3  ) 

ELSEIF  (ID.EQ.2)  THEN 
WX  =  WPER 
PLX  =  1.0 
ENDIF 

SEDQ(ISED)=(1.+SUSP)*SEDQ(ISED)*WX*PS 
TSEDQ(ISED)=TSEDQ(ISED)+SEDQ(ISED)*PLX*DTIM 
SUMPS=SUMPS+SEDQ( ISED) 

111  CONTINUE 

C  CALCULATE  THE  TOTAL  SEDIMENT  TRANSPORT  CAPACITY 

C  FOR  ALL  SIZES. 

SSEDQ=SSEDQ+SUMPS*DTIM*PLX 
110  CONTINUE 
121  CONTINUE 

IF(ID.LE.I)  THEN 

C  CALCULATE  THE  DETACHED  SEDIMENT  FROM  OVERLAND  FLOW. 

DCOH=DCOH*DTIM*PLNGTH(3 )*PLNGTH( I ) 

DA=DOF*( SSEDQ  -  TSDLAT) 

ELSEIF  (ID.EQ.2)  THEN 

C  CALCULATE  THE  DETACHED  SEDIMENT  FROM  CHANNEL  FLOW. 

DCOH=DCOH"'DTIM*PLNGTH(  3 ) 

DA=CHDOF* ( S  SEDQ-TSDLAT ) 

ENDIF 

IF  (DA.LT.0.0)  DA=0.0 

I F ( PLAS I ( I ) . GT . 0 . 0 , AND . DA . GT . DCOH )  DA=DCOH 
IF(ID.LE.l)  DA  =  TSDLAT  +  DA 
C 

C  REDISTRIBUTE  THE  AVAILABLE  SEDIMENT  IN  PROPORTION  TO  THE 

C  PARENT  MATERIAL  THEN  COMPUTE  THE  YIELD  FOR  EACH  SIZE. 
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TYIELD(I)  =  0.0 
DO  220  ISED-1,NSED 
IF(ID.LE.I)  THEN 
SEDX  =0.0 

ELSEIF  (ID.EQ.2)  THEN 
SEDX  =  SEDLAT(ISED) 

ENDIF 

SEDAV=SEDX  +  PSI(ISED)*DA 
YIELD(ISED, I )  =  TSEDQ(ISED) 

IF(TSEDQ( ISED) ,GT. SEDAV )  YIELD( ISED, I ) 
TYIELD(I)  =  TYIELD(I)  +  YIELD(ISED,I) 
220  CONTINUE 

RETURN 
END 


=  SEDAV 
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c 

c 


SUBROUTINE  POWER  ( Z , A , XJ I , XJ2 , CONV ) 


C 

C  THIS  SUBROUTINE  EVALUATES  THE  J1  AND  J2  INTEGRALS  IN  THE 

C  EINSTEIN  APPROXIMATION. 

C  XJ1  AND  XJ2  EQUAL  THE  J1  AND  J2  INTEGRAL  RESPECTIVLY 

C  N=ORDER  OF  APPROXIMATION  +1 

C  CONV=  CONVERGENCE  CRITERION 

N=  1 

XJ1=0. 

XJ2=0 . 

ALC=ALOG( A) 

C=  1 . 

D=-Z 
E  =  D+ 1 . 

FN  =  1 . 

AEX=A*--E 
GO  TO  102 

101  N=N+1 

C=C*D/ fn 

D=E 

E=D+1. 

FN=FLOAT(N) 

AEX=A**E 

102  IF  (ABS(E).LE. 0.001)  GO  TO  103 

XJ1=XJ1+C*( 1 .-AEX)/E 
X J  2 =X J  2 +C* ( ( AEX- 1 . ) / E**2  • -  AEX*ALG / E ) 

GO  TO  104 

103  XJ1=XJ1-C*ALG 

XJ2=XJ2-0 . 5*C*ALG**2 

104  IF  (N.EQ.l)  GO  TO  103 

CJ1=ABS( 1.-FJ1/XJ1) 

CJ2=ABS( 1.-FJ2/XJ2) 

IF  (CJ l.LE. CONV. AND. CJ2.LE. CONV)  .RETURN 

105  FJ1=XJ1 

FJ2=XJ2 
GO  TO  101 
END 
C 
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************************************************************** 


SUBROUTINE  OUT( I S EG , YI ELD , Q , NUM , QT , TAREA , TVINTR ,  NSED , DTIM , 

+  DMBI ,TYIELD,TRVOL, ID) 

C 

C  THIS  SUBROUTINE  CONTROLS  OUTPUT  TO  THE  SCREEN  AND  TO 
C  THE  DISK  FILE  SCREENl.OUT. 

C 

COMMON / VAR 1/  TITL 

DIMENSION  YIELD( 10) ,Q(201 ) ,DMBI ( 10) ,TITL( 3) 

C 

WRITE(6,99) 

WRITE(*,99) 

99  FORMAT (//////////) 

WRITE(6 , 100 ) I SEC, (TITL( I), 1=1, 3) 

WRITE (*, 100) I SEC, (TITL(I),I=1,3) 

100  F0RMAT(1X, 'THE  RESULTS  FOR  UNIT  NUMBER:  ' ,13, ' ( ' ,3AA, ' )' ) 
WRITE(6, 101 )TAREA,TRVOL, TVINTR, QT,TYIELD 

WRITE ( * , 1 0 1 )TAREA , TRVOL , TVINTR , QT , TYI ELD 

101  F0RMAT(1X,'  THE  TOTAL  AREA  IN  ACRES:  ' ,G12.5/, 

+  '  THE  TOTAL  RAINFALL  IN  ACRE  FEET:  ',G12.5,/, 

+  '  THE  TOTAL  INTERCEPTED  VOLUME  IN  ACRE  FEET:  ' ,G12.5 ,/, 

+'  THE  TOTAL  DISCHARGE  IN  ACRE  FEET:  ' ,G12.5,/, 

+  '  THE  TOTAL  SEDIMENT  YIELD  IN  POUNDS:  ',G12.5 ///) 

IF  (ID.LE.1)  THEN 
WRITE(6, 110) 

WRITE(*,110) 

110  FORMATC  (NOTE:  THE  FOLLOWING  VALUES  ARE  PER  FOOT  OF  WIDTH', 
+  '  FOR  A  PLANE  UNIT.)' ) 

ENDIF 

WRITE( 6 , 102 ) 

WRITE(*, 102) 

102  F0RMAT(/1X, 'THE  SEDIMENT  YIELD  BY  SIZES'/ 

+  '  SIZE  IN  MM  POUNDS' ) 

DO  AAA  ISED=1 ,NSED 

DMBIOUT  =  DMBI(ISED)*30A.8 
WRITE(6,103)  DMBIOUT, YIELD( I SED) 

WRITE(*, 103)  DMBIOUT, YIELD(ISED) 

AAA  CONTINUE 

103  F0RMAT(1X,G12.5,AX,G12.5) 

WRITE(6, 10A) 

WRITE(*,10A) 

10A  FORMATC ////IX, '  FINAL  HYDROGRAPH' //9X, 'TIME' 

-r  ,9X, '  DISCHARGE'  /8X, '  (MIN.  )  '  , 9X , '  (CFS )  '  ) 

DO  106  1=2, NUM 

TIME=(I-l)*DTIM/60. 

WRITE(6 , 105 )  TIME,Q(I ) 

WRITE(*, 105 )  TIME,Q(I ) 

105  FORMAT  ( 7X,F8 .2 ,6X,F10 . A ) 

106  CONTINUE 
RETURN 
END 


<cpyr>  5:  cat  code2 
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.MT  6 
.PO  12 
.OP 

PROGRAM  MUI.TSED3 
C 

C  THIS  IS  THE  CHANNEL  PROCRAM  COMPANION  TO  THE  MSED1  UPLAND 
C  WATERSHED  PROGRAM  MSED1 .  THIS  PROGRAM  MAKES  A  CALL  TO 

C  A  SUBROUTINE  CALLED  MSED2  WHICH  REORDERS  THE  OUTPUT  FROM  THE 

C  MSED1  PROGRAM  SO  THAT  THEY  CAN  BE  USED  IN  THIS  PROGRAM. 

C  THIS  PROGRAM  WAS  ORIGINALLY  DEVELOPED  AT  COLORADO  STATE  UNI- 

C  VERSITY  IN  THE  1970 ' S  AND  HAS  SUBSEQUENTLY  UNDERGONE  MODIFICA 
C  TIONS  BY  RESEARCHERS  AT  NEW  MEXICO  STATE  UNIVERSITY  AND  THE 
C  UNIVERSITY  OF  ILLINOIS. 

C  LAST  UPDATE:  JUNE  1987 

C 

C  THIS  PROCRAM  ROUTES  WATER  AND  SEDIMENT  FROM  UPLAND  WATERSHEDS 

C  AND  ADJACENT  CHANNEL  SIDE  SLOPES  THROUCH  A  CHANNEL  NETWORK. 

A  FINITE  DIFFERENCE  SOLUTION  OF  THE  METHOD  OF  CHARACTERISTICS 
FOR  A  KINEMATIC  WAVE  ASSUMPTION  IS  USED. 

C 

COMMON  /VAR1 /  NSEC , NWS , NPL, NDX, IMAX ,MM ,QOUT( 200 ) ,QDUM( 35 ) 
COMMON  /VAR2/  SEC( 35 ) , SLEN( 35 ) , SLOP1 ( 35 ) , DMB( 35 , 10 ) , 

1  PB0(35,10),DNB(35,10),WET  K( 35 ) , POROS( 35 ) , 

2  SI(35),SW(35), SUC(35 ) , VFALL(35 ,10) ,SAREA(35 ) 
COMMON  /VAR6/  NUM, SIN( 35 , 5 ) ,TSEDO( 35 , 10 ) ,QTOUT( 35 ) ,TIME(200 ) , 

1  TDDEP( 35 , 5 ) , ITYPE(35 ) 

COMMON  /VAR17/  TSEDEP( 35 ) , VCAP( 35 ) , VITL(35 ) , VT( 35 ) , 

1  CBUPIN(35,10),QIN(35),IFULL(35) 

DIMENSION  QSED(10),PPPM(10),SCOUR(35) 

CHARACTERS  FNAME9 
C 

C  CALL  MSED2  TO  SET  UP  THE  DATA  FILES  FOR  THIS  PROGRAM 
CALL  MSED2 

C  GET  SOME  INFORMATION  FOR  THIS  PROGRAM 
WRITE(*,2000) 

WRITE(*f 2001 ) 

READ(*,2004)  FNAME9 

2000  FORMAT(30(/),32X, 'WELCOME  TO  MSED3',/, 

121X, ' PC  VERSION  V87.06  BY  TIM  J.  WARD,  NMSU',15(/)) 

2001  FORMAT(10X, 'WHAT  IS  THE  NAME  OF  THE  MSED3  DATA  FILE?  '\) 

2004  FORMAT(A) 

C 

WRITE(*,2010)  FNAME9 

2010  FORMAT(/,20X, 'OPENING  FILE  ',A20,'  AS  FILE  9',/, 

120X, 'OPENING  FILE  SCREEN3 .OUT  FOR  SCREEN  OUTPUT',/, 

220X, 'OPENING  FILE  MSED3.DAT  AS  FILE  3',/, 

320X, 'OPENING  FILE  MSED4.DAT  AS  FILE  4',/, 

420X, 'OPENING  FILE  MSEDll.OUT  AS  FILE  11',/, 

520X, 'OPENING  FILE  MSED13.DAT  AS  FILE  13’,/, 

620X, 'OPENING  FILE  MSED14.DAT  AS  FILE  14',/) 

C 

OPEN (6 , FILE= ' SCREEN3 .OUT' ,STATUS='NEW' ) 

OPEN(ll,FILE=' MSEDll.OUT' , STATUS=' NEW ' ) 
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c 

c 

c 


900 

C 


C 


GO  GET  THE  DATA 
CALL  DATA1(FNAME9) 
ADJUST  THE  VISCOSITY 
CALL  TEMP 

INITIALIZE  THE  UNITS 
CALL  INITL 
WRITE(11,900) 
FORMAT(IX,'  ISEG 
*  /' 

CO  DO  THE  WORK 
CALL  ROUTE 
REWIND  6 
REWIND  11 
CLOSE(6 ) 

CLOSE ( 11 ) 

STOP  ' 

END 


AND  VARIABLES 


TIME  DISCHARGE 

(MINUTES)  (CFS) 


MSED3  IS  FINISHED’ 


SEDIMENT' , 
(PPM)' ) 
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c 

SUBROUTINE  MSED2 
C 

$ LARGE 

C  THIS  SUBROUTINE  REFORMATS  THE  FILES  CREATED  BY 
C  MSED1  FOR  USE  BY  MSED3 
C 

DIMENSION  Q ARRAY (200,70), SDARAY( 10 , 70 ) ,QT( 70 ) 
DATA  I READ , IWRIT1 , IWRIT2/7,3,4/ 

C 

WRITE/*, 2000) 

2000  FORMAT ( 15(/),3IX, 'MSED2  IS  RUNNING  *  ,  I5( / ) ) 
OPEN ( 3, F I LE=' MSED3.DAT* , STATUS= ' NEW ' ) 
OPEN(4,FILE=,MSED4. DAT 1 , STATUS= ' NEW ' ) 

OPEN( 7, FI LE=' MSED7.DAT' ) 

10  CONTINUE 

WRITE(*,2001)  IREAD 

2001  FORMATC  READINC  FILE  =  ',12/) 

WRITE(* ,2002 )  rWRITl , IWRIT2 

2002  FORMATC  WRITING  FILES  =  '  ,  I2,2X, 12///) 

READ(IREAD, 100)  NU,NUM1 

100  FORMAT/ 2 1 10) 

NUM=NUM1- 1 

DO  105  1=1, NU 

READ/ IREAD, 101)NS£D 

101  FORMAT/ 10X, 1 10) 

READ/ IREAD, 102 )(SDARAY/J, I ),J=1,NSED) 

102  FORMAT/ 4G 15 .6) 

WRITE/IWRIT2,102)/ SDARAY /J,l),J=l,NSED) 
READ/ IREAD ,102)QT(l) 

READ/ IREAD, 102 )/QAP RAY/ J, I ),J=1,NUM) 

105  CONTINUE 

WRITE/IWRIT1,102)/QT/I),I=1,NU) 

DO  110  J=1,NUM 

WRITE/lWRITl,102 ) (QARRAY (J,l),I=l,NU) 

110  CONTINUE 

IF/ IREAD. EQ. 8)  CO  TO  200 

IREAD=8 

IWRIT1= 14 

IWRIT2=13 

REWIND  3 

REWIND  4 

REWIND  7 

CLOSE/ 3) 

CLOSE (4 ) 

CLOSE/ 7) 

OPEN/8, FILE='MSED8. DAT' ) 

OPEN/ 1 3, FI LE=' MSED13.DAT' ,STATUS='NEW' ) 
0PEN(14,FILE='MSED14. DAT' , STATUS= ' NEW ' ) 

GO  TO  10 
200  CONTINUE 
REWIND  8 
REWIND  13 
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REWIND  14 
CL0SE(8) 
CLOSE (13) 
CLOSE (14) 
END 


c 

c 


SUBROUTINE  DATAI ( FNAME9 ) 


VoWoV 


c 

C  THIS  IS  WHERE  THE  DATA  ARE  READ  INTO  THE  PROGRAM 
C 

COMMON  / VAR1 /  NSEG . NWS ,NPL, NDX , IMAX ,MM , QOUT( 200 ) ,QDUM( 35 ) 
COMMON  /VAR2/  SEG( 35 ) , SLEN( 35 ) , SLOP1 ( 35 ) ,DMB( 35 , 10 ) , 

1  PB0(35,10), DNB(35 ,10) ,  WET  K(35 ) , POROS! 35 ) , 

2  SI(35),SW(35),SUC(35) , VFALL( 35,10), SAREA( 35 ) 
COMMON  /VAR3/  ISEC( 35 ) , IWS( 35 , 3 ) , IPL! 35 , 2 ) , IUP( 35 , 3 ) 

COMMON  /VAR4/  ACB( 35 ) , BEX( 35 ) , ADF( 35 ) , DELTS( 35 ) 

COMMON  /VAR5 /  VISCO! 35 ) ,T(35 ) 

COMMON  /VAR6/  NUM , S IN( 35 , 5 ) , TSEDO( 35 , 10 ) ,QTOUT( 35 ) , TIME( 200 ) , 

1  TDDEP( 35,5),ITYPE(35) 

COMMON  /VAR8/  A( 35 , 10 ) ,Q( 100) ,QP! 35 , 10) ,PB1( 35 , 10) ,GBO( 35 , 10 ) 
COMMON  / VAR9/  PORB( 35 ) , ZSUM , ASUM , EPS , DT , DX( 50 ) , DMAX( 35 ) 

COMMON  /VAR10 /CPR, EPR , ALAT ,QLAT , Al ( 35),B1(35),A2(35),B2(35) 
COMMON  /VAR  1  5  /  QWS  (  35  )  ,QPI.(  70  )  ,  QTWS  (  35  )  ,  QTPL(  70  )  ,  SEDPL(  70 , 10  )  , 
1  SEDWS( 35 , 10 ) ,CBOWS( 35 , 10) ,CBOPL( 70,10) 

COMMON  /VAR16/  XN(35) 

COMMON  /VAR17/  TSEDEP( 35 ) , VCAP( 35 ) , VITL( 35 ) , VT( 35 ) , 

1  CBUPIN(35,10),QIN(35), I  FULL! 35 ) 

COMMON  /VAR18/  PLAS I ( 35 ) , COHM( 35 ) 

DIMENSION  TITLE(20),PF(35, 1 1 ) , D< 35 , 11) 

CHARACTER'" 20  FNAME9 

n 

OPEN( 9 , FILE=FNAME9 , STATUS= ' OLD ' ) 

C 

C  INPUT  AND  OUTPUT  TITLE 

READ (9 ,110)(TITLE!lI),II=l,18) 

110  FORMAT( 20A4 ) 

WRITE ( 1 1 , 1 1 5 , ERR= 1 ) (TITLE( I I ) , I I =1 , 18 ) 

115  FORMLAT( '  TITLE:  ’,20A4) 

READ(9 , 150 )  DTIM, FTIM 
NUM=IFIX( FTIM/DTIM+ .00001 ) 

C  READ  IN  THE  CHANNEL  AND  RESERVOIR  INDICES 
READ! 9, 140)  NPL,NWS,NRES,NCH,NSED 
140  FORMAT! 81 10) 

C  SET  THE  MAXIMUM  NUMBER  OF  ITERATIONS  IN  WROUT 
C  AND  THE  NUMBER  OF  SECTIONS  THAT  THE  CHANNEL  IS  TO 
C  DIVIDED  INTO  FOR  ROUTING. 

IMAX=20 

NDX=5 

MM=NSED-1 

NSEC=NCH+NRES 

EPS=.01 

DT=DTIM 

DO  142  1=1, NSEG 

READ! 9, 145)  I SEG( I ) , IWS( 1 , 1 ) , IWS( 1 ,2 ) , IWS! 1 , 3 ) , IPL( 1 , 1 ) , 

1  IPL(I,2),IUP(I,1),IUP(I,2),IUP(I,3) 

142  CONTINUE 
145  FORMAT  (1015) 

DO  200  1=1, NSEG 
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READ(9, 160)  ITYPE(I) 

160  FORMAT(//I10) 

IF(ITYPE(I).EQ.1)  THEN 
C  THIS  IS  FOR  THE  CHANNEL  UNITS 

READ(9,150)  SLEN(I),SL0P1(I),WET  K(I),P0R0S(I), 

+  SI(I),SW(l), SUC(I ) , PLASlC I ) 

READ(9, 150)  XN(I),Td)  »COHM(I ) 

READ(9 , 150 )  A1(I),B1(I),A2(I),B2(I),ADF(I) 

C  SET  THE  TWO  MPM  PARAMETERS  AND  SHIELDS'  COEFFICIENT  TO  CONSTANTS 
AGB(I)  =  0.056 
BEX(I)  =  1.5 
DELTS(I)  =  0.047 
DO  180  J=1 ,NSED 

READ(9, 150)  D(I,J)fPF(I,J) 

180  CONTINUE 

ELSE 

C  THIS  IS  FOR  THE  RESERVOIR  UNITS 

READ (9 ,150)  VCAP{ I ) , VITL( I ) , SAREA( I ) , POROS( I ) 

VCAP( I )=VCAP(I )*43560 . 

VITL(I)=VITL(I)*43560. 

SAREA(I )=SAREA(I )*43560 . 

DO  187  J=1,NSED 

READ(9 , 150 )  D(I,J) 

187  CONTINUE 

ENDIF 

150  FORMAT(8F10.0) 

DO  185  J=1,MM 

IF(ITYPEd).EQ.l)  PB0(I ,  J  )*ABS(PF(I ,  J+l  )-PF(  I ,  J)  ) 
DNB(I,J)=(D(I,J+1)*D(I,J))**.5 
185  CONTINUE 
200  CONTINUE 
CLOSE(9) 

IF(NPL.EQ.O)  GO  TO  213 
OPEN ( 3 , FILE= ' MSED3 . DAT ' ) 

OPEN ( 4 , FILE= ' MSED4 . DAT ' ) 

DO  211  1=1, NPL 

READ(4,300)(SEDPL(I , JJ) , JJ=1 ,MM) 

211  CONTINUE 

READ (3 , 300 ) (QTPL( I ) , 1=1 ,NPL ) 

300  F0RMAT(4G15.6) 

CLOSE(4) 

213  CONTINUE 
IF(NWS.EQ.O)  GO  TO  214 
OPEN( 13, FILE= ' MSED1 3 . DAT ' ) 

OPEN( 14, FILE=' MSED14.DAT' ) 

DO  212  1=1, NWS 

READ( 13 ,300)(SEDWS( I , JJ) , JJ=1 ,MM) 

212  CONTINUE 

READ( 14,300) (QTWS(I ) , 1=1 , NWS) 

CLOSE( 13) 

214  CONTINUE 
RETURN 

1  STOP  '  **********  ERROR  ON  WRITE  TITLE  ********** ' 

END 


C 
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SUBROUTINE  TEMP 
C 

C  -  THIS  SUBROUTINE  CORRECTS  THE  VISCOSITY  AND 

C  -  HYDRAULIC  CONDUCTIVITY  FOR  TEMPERATURE  VARIATIONS 

C  —  FROM  THE  ASSUMED  TEMPERATURE  OF  68  DEGREES  (F). 

C  -  PARAMETER  DEFINTIONS . 

C  T  TEMPERATURE  IN  DEGREES  F. 

C  VISCO  =  KINEMATIC  VISCOSITY  (FT**2/SEC) 

C  I  PLANE  =  NUMBER  OF  PLANES. 

C 

COMMON  / VAR1 /  NSEG , NWS ,NPL , NDX , IMAX ,MM ,QOUT( 200 ) ,QDUM( 35 ) 

COMMON  /VAR2/  SEG( 35 ) , SLEN( 35 ) , SLOP1 ( 35 ) , DMB( 35 , 10 ) , 

1  PB0(35,10),DNB(35,10),WET  K(35 ) ,POROS( 35 ) , 

2  SI(35),SW(35),SUC(35),VFALL(35,10),SAREA(35) 
COMMON  /VAR5/  VISCO( 35 ) ,T( 35 ) 

DIMENSION  TE(10),V(10) 

DATA  TE/ 32.  ,40.  ,50.  ,60. ,68.  ,80.  ,90.  ,100.  ,120.,  140./, 

1  V/ 1.93, 1.66, 1.41, 1.22, 1.09, 0.930, 0.826, 0.739, 0.609, 0.5 14/ 

C 

C  —  CALCULATE  NEW  VISCOSITY  BY  INTERPOLATION. 

DO  110  K=1,NSEC 
DO  100  1=1,10 

IF(TE(I).LT.T(K))  GO  TO  100 

FAC1=(T(K)-TE(I-1))/(TE(I)-TE(I-1>) 

VISC0(K)=V(I-1)+FAC1*(V(I)-V(I-1)) 

GO  TO  110 

100  CONTINUE 
110  CONTINUE 

C  —  ADJUST  THE  HYDRAULIC  CONDUCTIVITY. 

DO  101  K=1,NSEG 
10  FAC2=VISCO(K)/ 1.09 

WET  K(K)=WET  K(K)/FAC2 
VI SCO(K)=VISCO(K)*. 00001 

101  CONTINUE 
RETURN 
END 

C 
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SUBROUTINE  INITL 


THIS  SUBROUTINE  IS  USED  TO  INITIALIZE  THE  VARIABLES  AND  SET  UP 
THE  CORRECT  UNITS  FOR  FURTHER  COMPUTATIONS 


COMMON 

COMMON 

1 

2 

COMMON 

COMMON 

COMMON 

COMMON 

1 

COMMON 

1 

COMMON 

COMMON 

COMMON 

1 


/VARl/  NSEG,NWS,NPL,NDX,IMAX,MM,QOUT(200),QDUM(35) 

/VAR2/  SEG(35),SLEN(35),SLOP1(35),DMB(35,10), 

PB0(35 , 10) , DNB(35 , 10) , WET  K(35) ,POROS(35), 
SI(35),SW(35),SUC(35),VFALL(35,10),SAREA(35) 

/VAR3/  ISEG(35) ,IWS(35 ,3) , IPL(35 , 2) ,IUP(35 , 3 ) 

/VAR4/  AGB(35 ) ,BEX( 35 ) ,ADF( 35 ) , DELTS(35 ) 

/VAR5/  VISCO(35 ) ,T( 35  ) 

/VAR6/  NUM, SIN(35 , 5 ) ,TSEDO( 35 , 10) ,QTOUT( 35 ), TIME (200) , 
TDDEP(35 ,  5 ) , ITYPE( 35 ) 

/VAR7/  EGB( 35 ) ,GBC( 35 ) ,ZBL( 35 ) ,GBUP( 35 ) ,GBLAT(35 ) ,BLAT(35 ) , 
BTEM( 35), SB (35, 5, 10), RB (35, 5, 10) f SLOP(35 , 5 ) ,DDEP(5 ) 
/VAR8/  A(35,10),Q(100),QP(35,IO),PB1(35,10),GBO(35,!0) 
/VAR9/  PORB(35 ) ,ZSUM, ASUM, EPS , DT,DX(50) ,DMAX(35) 

/VAR17/  TSEDEP( 35 ) , VCAP( 35 ) ,VITL( 35 ) , VT(35 ) , 
GBUPIN(35,10),QIN(35),IFULL(35) 


INITIALIZE  ENTIRE  WATERSHED 

ITC0M1=NUM 

DO  100  M=1 ,MM 

EGB(M)=0. 

GBC(M)=0. 

ZBL(M)=0. 

BTEM(M)=0. 

100  CONTINUE 

DO  300  N=1 ,NSEG 
DX(N)=SLEN(N)/5. 

QIN(N)=0.0 

IFULL(N)=0 

TSEDEP(N)=0.0 

IF(ITYPE(N).EQ.2)  VT(N)=VITL(N) 

Q(N)=0. 

IF(ITYPE(N).EQ.1)WET  K(N)=WET  K(N)/43200. 
QTOUT(N)=0. 

QDUM(N)=0. 

DO  300  J=1,NDX 
GBUPIN(N,M)=0 .0 
A(N, J)=0. 

QP(N,J)=0. 

SIN(N, J)=0 . 

IF( ITYPE(N) .EQ. 1 )SLOP(N, J)=SL0P1(N) 
TDDEP(N, J)=0. 

DO  300  M=1 ,MM 
GBO(N,M)=0. 

TSEDO(N,M)=0. 

SB(N,J,M)=0. 

RB(N,J,M)=0. 

300  CONTINUE 
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C  CONVERT  SEDIMENT  SIZES,  FIND  FALL  VELOCITIES  AND  MAXIMUM 

C  PARTICLE  SIZE,  AND  CALCULATE  FRACTION  OF  THE  SOIL  THAT  IS  SOLID. 

DO  106  N=1 , NSEG 
DO  104  M=1,MM 
DMB(N,M)=DNB(N,M)/304.8 
IF  (DMB(N,M)  .LE.  0.0002)  THEN 

VFALL( N , M ) =2 . 95 1 7*DMB ( N , M)**2/ VI SCO ( N ) 

ELSE 

VFALL(N ,M)=( SQRT( 36 .064'vDMB(N ,M)**3+36 ,*VISCO(N)**2 )-6 .*VISCO(N) ) 
1  /DMB(N,M) 

ENDIF 

104  CONTINUE 

DMAX (N)=DMB(N,MM) 

PORB( N )=1 . -POROS ( N ) 

106  CONTINUE 
RETURN 
END 
C 
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SUBROUTINE  ROUTE 

THIS  IS  THE  SUBROUTINE  THAT  DOES  ALL  OF  THE  WORK. 

COMMON  /VARl/  NSEG,NWS,NPL,NDX, IMAX,MM,QOUT(200) ,QDUM(35 ) 

COMMON  /VAR2/  SEG(35),SLEN(35),SLOP1(35),DMB(35,10), 

1  PB0(35 , 10) ,DNB(35 , 10) ,WET  K(35) ,POROS( 35 ) , 

2  SI(35),SW(35), SUC( 35 ) ,VFALL(35 ,10) , SAREA(35 ) 

COMMON  /VAR3/  ISEG(35 ) ,IWS(35 , 3 ) , IPL(35 ,2) ,IUP(35 , 3 ) 

COMMON  /VAR4/  AGB(35) ,BEX( 35 ) , ADF( 35 ) ,DELTS( 35 ) 

COMMON  /VAR5/  VISCO( 35 ) ,T( 35 ) 

COMMON  /VAR6/  NUM,SIN( 35 ,5 ) ,TSEDO(35 , 10 ) ,QTOUT(35 ) ,TIME(200) , 

1  TDDEP(35 , 5 ) , ITYPE( 35 ) 

COMMON  /VAR7 /  EGB( 35 ) ,GBC( 35 ) , ZBL( 35 ) , GBUP( 35 ) ,GBLAT( 35 ) , BLAT( 35 ) , 
1  BTEM(35 ), SB (35, 5, 10), RB (35, 5, 10) , SLOP (35 , 5 ) ,DDEP( 5 ) 

COMMON  /VAR8/  A(35 , 10 ) ,Q( 100 ) ,QP( 35 , 10 ) , PB1(35 , 10) ,GBO(35 , 10) 
COMMON  /VAR9/  PORB(35) ,ZSUM, ASUM, EPS ,DT,DX(50) ,DMAX(35 ) 

COMMON  /VAR10/  CPR,EPR, ALAT,QLAT ,Al( 35 ) ,B1( 35 ) , A2(35 ) ,B2(35 ) 

COMMON  /VAR16/  XN(35) 

DIMENSION  QSED(10),PPPM(10) 

C 

DTS=DT*60 . 

DTN=DTS*FLOAT ( NDX ) 

ITCOMl=NUM 

DO  100  IT=1 , ITCOM1 

C  IF  APPROPRIATE,  READ  FILES  FROM  MSED1 
IF(NWS.GE. 1 .OR.NPL.GE. 1 )CALL  DATA2 
TIME(IT)=FLOAT(IT)*DT 
C  COMPUTE  AT  TIME  IT  (T+DT) 

DO  200  1=1 ,NSEG 
K=ISEG( I ) 

WRITE(*,2001)  I, IT 

2001  FORMAT(2X, 'SEGMENT=  ' ,15 ,5X, 'TIME  STEP=  ’,15) 

IF(ITYPE(K) . EQ.2)  THEN 
CALL  RES(K,QUP,DTS) 

Q(K)=QUP 
GO  TO  311 
END  IF 
CPR=A1 (K) 

EPR=B1(K) 

DTX=DTN/ SLEN(K) 

C  DETERMINE  THE  UPSTREAM  AND  LATERAL  CHANNEL  INFLOWS 
CALL  UPLAT  (K,DTS,SLEN(K) ,QUP) 

DO  300  J=1 ,NDX 
BPER=A2  (  K ) 

ZSUM=0. 

SLP=SLOP(K, J) 

DO  310  M=1 ,MM 
ZBL(M)=SB(K,J,M) 

ZSUM=ZSUM+ZBL(M) 

BTEM (M)=RB(K,J,M) 

310  CONTINUE 
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QCONC=QUP+QLAT*SLEN( K) /NDX 
Q I  .AT P  =  QI.AT 

C  IF  SOME  WATER  IS  PRESENT,  TRY  TO  INFILTRATE  IT 

IF(A(K,J).GT. 1.0E-04)  CALL  CHINL(K , SIN(K , J ) ,QUP , DTS , A(K , J ) ,QLATP) 

ALAT=QLATP'-DTS 

ASUM=ALAT+A( K , J )+DTX*QUP 

C  DETERMINE  THE  COEFFICIENT  AND  THE  EXPONENT  IN  A-Q  RELATION 
CALL  RES  I ST( XN( K ) , A1 ( K ) , B I ( K ) , SLP , ALP , BET ) 

IF(ASUM. LE.  1.0E-05)  GO  TO  320 
C  IF  WE  HAVE  SOME  WATER  LEFT,  ROUTE  IT  IN  THE  CHANNEL 

CALL  WROUT  ( K , J , DTX , ALP , BET , QUP ,QE ) 

IF(QE.LE.O.OI)  GO  TO  320 
C  REDISTRIBUTE  THE  SEDIMENT  IN  THE  CHANNEL 
CALL  PERTG(K) 

C  IF  WE  STILL  HAVE  ENOUGH  FLOW  TO  USE,  TRANSPORT  SOME  SEDIMENT 
C  AND  ROUTE  IT  THROUGH  THE  CHANNEL 

CALL  TRANSP  ( K , J , BPER , SLP , QE , AREA, ALP , BET, DC , DCOH ) 

CALL  SROUT  ( K , J , DTX , BPER , AREA ,QE , DC , QCONC , DCOH ) 

GO  TO  350 

C  IF  NO  MORE  (NOT  MUCH)  WATER,  TAKE  CARE  OF  THE  ROUTING  IN  CEASE 
320  CALL  CEASE  (K , J , DTX , QUP, QE, AREA) 

350  A(K,J)=ALP*QE**BET 
QP(K, J )=QE 
QUP=QE 
SUMV=0 . 

DO  390  M=L,MM 

GBUP( M) =RB( K, J ,M)*QE 

SB(K,J,M)=ZBL(M)+EGB(M)/BPER 

GBO(K,M)=GBUP(M) 

sumv=sumv+egb(m) 

390  CONTINUE 

DDEP( J)=SUMV/(PORB(K)*BPER) 

TDDEP(K,J)=TDDEP(K,J)+DDEP(J) 

300  CONTINUE 

Q(K)=QUP 

311  CONTINUE 

IF( K . EQ . NSEG )QOUT( IT )=QUP 
GBOSUM=0 . 0 
DO  433  M=1 ,MM 

GBOSUM=GBOSUM+GBO(K,M)*2 .65 
433  CONTINUE 
PPM=0 . 

DO  312  M= 1 , MM 

QSED ( M )  =GBO( K , M )* 1 6  5 . 4* DTS 
PPPM(M)=0.0 

IF(QUP.GT.O.O)PPPM(M)=GBO(K,M) *2 650000 . / (QUP+GBOSUM) 
PPM=PPM+PPPM(M) 

312  CONTINUE 

CALL  OUT( K ,QSED ,QUP , PPM , PPPM , IT , ITCOM1 , DTS ) 

200  CONTINUE 
100  CONTINUE 
RETURN 
END 
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************************************************************** 
SUBROUTINE  DATA2 

THIS  SUBROUTINE  IS  USED  TO  READ  DATA  FILES  CREATED  BY  MSEDI 

COMMON  /VAR1/  NSEG,NWS ,NPL,NDX, IMAX,MM,QOUT(200 ) ,QDUM(35 ) 

COMMON  /VAR15/  QWS(35 ) ,QPL( 70 ) ,QTWS( 35 ) ,QTPL( 70) , SEDPL( 70, 10) , 

1  SEDWS(35 , 10) ,GBOWS(35 , 10) ,GBOPL( 70, 10) 

IF(NWS.LE.O)  GO  TO  120 
READ(14,100)(QWS(I),I=1,NWS) 

100  FORMAT(4G15.6) 

DO  110  1=1, NWS 

DO  110  JJ=1 ,MM 

IF  (QTWS(I).LE.O.O)  THEN 
GBOWS(I , JJ)=0 .0 
ELSE 

GBOWS( I , JJ )=SEDWS( I , JJ )*QWS ( I )/ (QTWS( I )*43560.*165 .4) 
END  IF 

110  CONTINUE 

120  CONTINUE 

IF(NPL.LE.O)  RETURN 

READ( 3 , 100 ) (QPL( I ) , 1=1 , NPL ) 

DO  200  1=1, NPL 

DO  200  JJ=1 ,MM 

IF(QTPLd).LE.O.O)  THEN 
GBOPL(I , JJ)=0.0 
ELSE 

GBOPL(I , JJ)=SEDPL(I , JJ)*QPL(I )/(QTPL(I )*43560.*165 .4) 
END  IF 

200  CONTINUE 
RETURN 
END 
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c 

c 

c 

c 

c 

c 

c 

c 


c 


****************** 


SUBROUTINE  RES( K ,QUP , DTS ) 

THIS  SUBROUTINE  IS  USED  TO  MOVE  WATER  AND  SEDIMENT  THROUGH 
A  RESERVOIR  UNIT.  THE  RESERVOIR  WILL  FILL  THEN  SPILL  AT  THE 
SAME  RATE  THAT  WATER  ENTERS.  SEDIMENT  WILL  BE  TRAPPED  ACCORDING 
TO  TRPEF .  MORE  EFFORT  NEEDS  TO  BE  EXPENDED  ON  THIS  SUBROUTINE 

DIMENSION  CBUPAV ( 10 ) ,QINPRE( 35 ) ,CBUPRE( 35 , 10  ) 

COMMON  / VAR  1 /  NSEG , NWS , NPL , NDX , IMAX ,MM , QOUT( 200 ) ,QDUM( 35 ) 

COMMON  /VAR2 /  SEC( 35 ) , SLEN( 35 ) , SLOP  1 ( 35 ) , DMB( 35 , 10 ) , 

1  PBO(35, 10) , DNB( 35 , 10) , WET  K( 35 ) , POROS(35 ) , 

2  ST(35),SW(35),SUC(35) , VFALL( 35,10), SAREA( 35 ) 

COMMON  /VAR 7!  EGB( 35 ) , GBC( 35 ) , ZBL( 35 ) ,CBUP( 35 ) ,GBLAT( 35 ) , BLAT( 35 ) , 

1  BTEM( 35),SB(35,5,10),RB(35,5,10), SLOP( 35 ,5),DDEP(5) 

COMMON  /VARS/  A( 35 , 10 ) ,Q( 100 ) ,QP( 35 , 10 ) , PB1 ( 35 , 10 ) ,CBO( 35 , 10 ) 
COMMON  /VAR9/  PORB( 35 ) , ZSUM, ASUM, EPS , DT, DX( 50 ) , DMAX( 35 ) 

COMMON  VAR  17/  TSEDEP( 35 ) , VCAP( 35 ) , VITL( 35 ) , VT( 35 ) , 

1  GBUPIN(35,10),QIN(35),IFULL(35) 


QINPRE(K)=QIN0C) 

C  FIND  THE  UPSTREAM  INFLOW  TO  THE  RESERVOIR 

CALL  UPRES(K.QUP) 

QIN(K)=QUP 

QAV=(QINPRE ( K)+QIN(K))/2. 

DO  100  M= 1 , MM 

CBUPRE(K,M)=CBUPIN(K,M) 

gbupin(k,m)=gbup(m) 

cbupav(mMcbup(m)+gbupre(k,m))/2. 

100  CONTINUE 


VIN=QAV'"DTS 

VT(K)=VT(K)+VIN 

I  F( VT(K) . LT.VCAP(K) )  GO  TO  400 

IFULL(K)=IFULL(K)+1 

QUP=QAV 

I F ( IFULL(K) . GE. 2 )  CO  TO  200 
QUP=(VT(K)-VCAP(K) )/DTS 
200  CONTINUE 

VOVFL=1.5*QUP/SAREA(K) 

IF(VOVFL.LE.O.O)  GO  TO  400 
DO  300  M= 1 , MM 

TRPEF=VFALL(K,M)/VOVFL 
IF(TRPEF.GT. 1 .0 )TRPEF=1 .0 
GBO(K,M)  =  (  1  .-TRPEF)'-'-GBUPAV(M)--QUP/QAV 
SEDDEP=(GBUPAV(M)-GBO(K,M)  )*DTS 
TSEDEP(K)=TSEDEP(K)+SEDDEP/ (43560. *PORB(K)) 
300  CONTINUE 
RETURN 

400  CONTINUE 


QUP=0 . 0 
DO  500  M= 1 , MM 

GBO( K , M ) =0 . 0 

SEDDEP=( GBUPAV (M)-GBO( X, M )  )'VDTS 
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TSEDEP(K)=TSEDEP(K)+SEDDEP/( 43560. *PORB(K) ) 
500  CONTINUE 
RETURN 
END 
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SUBROUTINE  Ul’RES  (  K  ,QUP  ) 


c 

C  THIS  SUBROUT  INK  TS  USED  TO  FIND  THE  UPSTREAM 
O  INFLOWS  TO  A  RESERVOIR  UNIT. 

C 

COMMON  /VAR  1  /  NSEC , NWS  ,  NPl. ,  NDX ,  IMAX ,  MM ,QOUT(  200 )  ,Q[)UM(  35  ) 

COMMON  /VAR 3/  l SEC( 35  ) , I WS( 35 , 3 ) , IPL( 35 , 2 ) , IUP( 35 , 3 ) 

COMMON  /VAR  7/  ECB( 35 ) ,CBC( 35 ) , ZBL( 35 ) ,CBUP<  35  )  ,GBLAT(  35  ) ,  BLAT(  35  )  , 
1  BTEM( 35 ) , SB( 35 , 5 , 10 ) ,RB( 33 , 5 , 10 ) , SLOP(  35 , 5 ) , DDEP(  5 ) 

COMMON  /VAR8/  A( 35 , 10 ) ,Q( 100 )  ,QP( 35 , 10 ) , PB1 (35 , 10 ) ,CBO( 35 , LO ) 
COMMON  /VAR  15/  QWS( 35 ) ,QPL( 70 ) ,QTWS( 35 ) ,QTPL( 70) , SEDPL( 70 , 10) , 

1  SE,;WS  (35,10),  CBOWS(  35,10)  ,CBOPL(  70,10) 

COMMON  / VAR  17/  TSEI)EP(  35  )  ,  VCAP(  35  ) ,  VITL(  35  )  ,  VT(  35 ) , 

1  CBUPrN(35, 10),QIN(35),rFULL(35) 

C 

QUP=0 . 

DO  101  M=1 ,MM 
GBUP(M )=0 . 0 
101  CONTINUE 
200  FORMAT ( / F 1 0 . 0 ) 

210  FORMAT( 8F 10 . 0 ) 

DO  100  J=1 , 3 

IF(  IUP(K,  J)  .F.Q.O)  CO  TO  105 
JJ=1UP(K, J) 

QUP-QUP+Q(JJ) 

DO  100  M=1 , MM 

CBUP(M)=CBUP(M)+CBO(JJ,M) 

100  CONTINUE 
105  CONTINUE 

DO  110  J=1 , 3 

IF(IWS(K,J)-EQ.O)  CO  TO  115 
JJ=IWS (K , J ) 

QUP=QUP+QWS( JJ ) 

DO  110  M=1,MM 

CBUP(M)=CBUP(M)+CBOWS( JJ,M) 

110  CONTINUE 
115  CONTINUE 
RETURN 
END 
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SUBROUTINE  UPLAT  (K,DTS,SLEN,QUP) 

THIS  SUBROUTINE  DETERMINES  THE  UPSTREAM  AND  LATERAL  INFLOW 
C  TO  A  CHANNEL  UNIT. 

C 

COMMON  /VAR1/  NSEG,NWS,NPL,NDX,IMAX,MM,QOUT(200),QDUM(35) 

COMMON  /VAR3/  ISEG( 35 ) , IWS(35 , 3 ) ,IPL( 35 , 2) , IUP(35 ,3 ) 

COMMON  /VAR7/  EGB( 35 ) ,GBC(35 ) ,ZBL( 35 ) ,GBUP(35 ) ,GBLAT( 35 ) ,BLAT(35 ) , 
1  BTEM(35) ,SB(35 ,5 , 10 ) ,RB(35 , 5 , 10) ,SLOP(35 , 5 ) ,DDEP( 5 ) 

COMMON  /VAR8/  A(35 , 1") ,Q(100) ,QP( 35 , 10) ,PB1(35 , 10 ) ,GBO(35 , 10 ) 
COMMON  /VAR9/  PORB(  3‘> )  ,ZSUM,  ASUM,  EPS , DT,DX(50)  ,DMAX( 35  ) 

COMMON  /VAR10/  CPR,EPR,ALAT,QLAT,A1 (35 ) ,B1 (35 ) ,A2 (35 ) ,B2( 35 ) 

COMMON  /VAR15/  QWS(35 ) ,QPL( 70 ) ,QTWS(35 ) ,QTPL( 70) ,SEDPL( 70 , 10) , 

1  SEDWS(35 ,10) ,GBOWS(35 ,10) ,GBOPL(70, 10) 

C 

QUP=0 . 

QLAT=0 . 

ALAT=0. 

C  DETERMINE  THE  UPSTREAM  INFLOW  RATE 
DO  11  1=1, MM 
BLAT( I )=0 . 

GBUP(I)=0. 

GBLAT(I)=0. 

11  CONTINUE 

DO  109  J=l,  3 

IF  (IUP(K,J)  .EQ.  0)  GO  TO  109 
JJ=IUP(K,J) 

QUP=QUP+Q(JJ) 

DO  108  M=1,MM 
GBUP(M)=GBUP(M)+GBO(JJ,M) 

108  CONTINUE 

109  CONTINUE 

DO  110  J=l,3 

IF(IWS(K,J) . EQ.O)  GO  TO  ill 
JJ=IWS(K, J) 

QUP=QUP+QWS ( J J ) 

DO  110  M=1 ,MM 

GBUP(M)=GBUP(M)+GBOWS(JJ,M) 

110  CONTINUE 

111  CONTINUE 

C  DETERMINE  THE  LATERAL  INFLOW  RATE 
DO  113  J=l,  2 

IF(IPL(K,J).EQ.O)  GO  TO  113 
JJ=IPL(K, J) 

QLAT=QLAT+QPL(JJ) 

DO  112  1=1, MM 

GBLAT( I )=GBLAT( I )+GBOPL( JJ, I )*SLEN 
BLAT( I ) =GBLAT( I )*DTS 

112  CONTINUE 

113  CONTINUE 

200  FORMAT(/F10.0) 

210  FORMAT(8F10.0) 
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RETURN 

END 
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C  ***********************VHt********T^*************,*TSr**********’** 

c 

SUBROUTINE  CHINL(K,SIN, QUP , DTS , AREA , QLATP ) 

C 

C  -  THIS  SUBROUTINE  CALCULATES  CHANNEL  INFILTRATION. 

C  -  PARMETER  DEFINITIONS. 

C  SI  INITIAL  SATURATION  OF  CHANNEL  SOIL. 

C  SW  WETTED  SATURATION  OF  CHANNEL  SOIL. 

C  SUC  =  AVERAGE  SUCTION  OF  CHANNEL  SOIL  (INCHES). 

C  POROS  =  POROSITY  OF  SOIL. 

C  WET  K  =  HYDRAULIC  CONDUCTIVITY  OF  CHANNEL  SOIL  (IN/HR). 

C  SIN  =  ACCUMALATED  INFILTRATED  VOLUME  (FT**3). 

C  QUP  =  UPSTREAM  FLOW  (CFS). 

C  QLATP  =  LATERAL  INFLOW  (CFS/FT). 

C  DTIM  =  LENGTH  OF  TIME  INTERVAL  (SEC). 

C  AREA  =  PREVIOUS  CROSS  SECTIONAL  AREA  (FT**2). 

C 

COMMON  /VAR1 /  NSEG,NWS,NPL,NDX, IMAX ,MM,QOUT( 200) ,QDUM( 35 ) 
COMMON  /VAR2/  SEG( 35 ) , SLEN( 35 ) , SLOP1 ( 35 ) ,DMB( 35 , 10 ) , 

1  PBO(35, 10) ,DNB(35 ,10) ,WET  K(35 )  ,POROS(35 ) , 

2  SI(35),SW(35), SUC(35 ) ,VFALL( 35 , 10) ,SAREA( 35 ) 
COMMON  /VAR9/  PORB(35 ) ,ZSUM, ASUM,EPS ,DT,DX(50 ) ,DMAX(35 ) 

COMMON  /VAR10/  CPR,EPR,ALAT,QLAT,A1(35) ,B1(35) ,A2(35),B2(35) 

C 

C  —  GO  BACK  IF  CONDUCTIVITY  WAS  SET  TO  0.0  TO  AVOID 

C  -  AN  INFILTRATING  CHANNEL. 

IF(WET  K(K).LE.O.O)  RETURN 

C  -  CALCULATE  CHANNEL  TOP  WIDTH  AND  DEPTH  OF  FLOW. 

TW=A2(K)*AREA**B2(K) 

D=AREA/TW 

C  CALCULATE  POTENTIAL  INFILTRATION. 

AL=(SUC(K)/ 12 ,+D)*(SW(K)-SI(K) )*POROS(K) 

C1=2.*SIN-WET  K(K)*DTS 
C2=8.*WET  K(K)*DTS*( AL+SIN) 

SEGL=SLEN ( K ) / FLOAT ( NDX ) 

DELF=( ( -Cl+( Cl**2+C2 )** .  5 )  /2 . )*SEGL*TW 
IF(DELF.EQ.O.O)  RETURN 
C  CALCULATE  VOLUME  OF  AVAILABLE  WATER. 

DELV= ( QLAT*S  EGL+QUP ) *DTS 
C  —  COMPARE  POTENTIAL  WITH  AVAILABLE. 

IF(DELF.GE.DELV)GO  TO  30 
SIN  =  SIN  +  DELF / ( SEGL*TW ) 

QLATP=QLAT-DELF / ( SEGL*DTS ) 

RETURN 

C  -  THIS  IS  THE  CASE  WHERE  ALL  THE  AVAILABLE 

C  —  WATER  IS  INFILTRATED. 

30  SIN  =  SIN  +  DELV/(SEGL*TW) 

QLATP=0. 

QUP=0 . 

RETURN 

END 

C 

Q  *****«r***TSr****1ir*i!ni-**T!r*TS-************************************** 

C 
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SUBROUTINE  RESIST(XN,A1 ,B1 ,SLP,ALP,BET) 

C 

C  THIS  SUBROUTINE  CALCULATES  THE  PARAMETERS  A  AND  B 
C  IN  THE  EQUATION  AREA=A*Q**B . 

C 

C  —  THIS  IS  FOR  THE  MANNING'S  RESISTANCE. 

BET=(3.)/(5.-2.*Bl) 

ALP=( Al** (4 . / 3 . )*XN**2 . / ( 2.2 1*SLP) )**(3 . /( 10.-4. *B1 ) ) 
RETURN 
END 
C 
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SUBROUTINE  WROUT  (K, J,DTX,ALP,BET,QUP,QE) 

C 

C  THIS  SUBROUTINE  DOES  THE  KINEMATIC  WAVE  WATER  ROUTING. 

C  THE  TECHNIQUE  ROUTES  AREAS  (DEPTHS),  NOT  DISCHARGES. 

C 

COMMON  /VAR1/  NSEG,NWS,NPL,NDX,IMAX,MM,QOUT(200),QDUM(35) 
COMMON  /VAR8/  A(35 , 10) ,Q( 100) ,QP(35 , 10) ,PB1( 35 , 10) ,GBO(35 , 10) 
COMMON  /VAR9/  PORB(35),ZSUM,ASUM,EPS,DT,DX(50),DMAX(35) 

COMMON  /VAR10/  CPR,EPR,ALAT,QLAT,A1(35),B1(35),A2(35),B2(35) 

C 

C  SET  UP  THE  A-Q  RELATIONSHIP 
QPRX=QP(K,J) 

QAVE=0 . 5* (QUP+QPRX ) 

BEM=BET-1 . 

BEN=BEM-1. 

ALBET=ALP*BET 

ALBEM=ALP*BET*BEM 

DTXA=DTX+ALP 

ERROR=EPS*ASUM 

C  LINEAR  SCHEME  TO  FIND  THE  FIRST  APPROXIMATION 
ITER=0 

IF  (QAVE.GT. 1 . OE-5)  THEN 
DAQ=ALBET*QAVE**BEM 
QE=( ALAT+DTX*QUP+DAQ*QPRX ) / ( DTX+DAQ ) 

ELSE 

QE=ASUM/DTXA 

ENDIF 

C  NONLINEAR  SCHEME  TO  REFINE  THE  SOLUTION 

115  ITER=ITER+1 
AEST=DTX*QE+ALP*QE**BET 
ADEV=ASUM-AEST 

IF  (ABS(ADEV).LE. ERROR)  GO  TO  120 
IF  (ITER.LT. IMAX)  GO  TO  116 

STOP'  **********  STOPPED  IN  WROUT  ********** ' 

116  FDER=DTX+ALBET*QE**BEM 
SDER  >LBEM*QE**BEN 
BB=FDER/SDER 

SC=2 . *ADEV / SDER 

STEM=BB*BB+SC 

IF  ( STEM . GE . 0 . )  GO  TO  117 

QE=QE+ADEV/FDER 

GO  TO  115 

117  STEM=SQRT ( STEM ) 

IF  (ADEV.GT.O . )  GO  TO  119 

ETEM=BB+STEM 

QE=QE-ETEM 

IF  (QE.GT.O.)  GO  TO  115 

118  ETEM=0. 5*ETEM 
QE=QE+ETEM 

IF  (QE.GT.O.)  GO  TO  115 
GO  TO  118 

119  X1=QE-BB-STEM 
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X2=QE-BB+STEM 

AD 1 =ABS ( ASUM- DTX*X 1 -ALP*X 1 **BET ) 
AD2=ABS(ASUM-DTX*X2-ALP*X2**BET) 
QE=X1 

IF  (AD1.GT.AD2)  QE=X2 
GO  TO  115 

120  RETURN 
END 
C 
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SUBROUTINE  PERTG(K) 

C 

C  THIS  SUBROUTINE  IS  USED  TO  REDISTRIBUTE  THE  SEDIMENT 
C  IN  THE  STREAM  AS  IT  IS  DEPOSITED  OR  ERODED. 

C 

COMMON  /VARl/  NSEG,NWS,NPL,NDX,IMAX,MM,QOUT(200),QDUM(35) 

COMMON  /VAR2/  SEG(35 ) ,SLEN( 35) , SL0P1(35 ) ,DMB(35 , 10) , 

1  PB0(35 ,10) ,DNB( 35 , 10 ) ,WET  K(35 ) ,POROS( 35 ) , 

2  SI(35),SW(35),SUC(35),VFALL(35,10),SAREA(35) 

COMMON  /VAR7/  EGB(35 ) ,GBC(35 ) , ZBL(35 ) ,GBUP(35),GBLAT(35) ,BLAT(35) , 

1  BTEM(35),SB(35,5,10),RB(35,5,10),SLOP(35,5),DDEP(5) 

COMMON  /VAR8/  A(35 , 10) ,Q( 100) ,QP( 35 , 10) ,PB1( 35 , 10) ,GBO(35 , 10 ) 
COMMON  /VAR9/  PORB(35 ) ,ZSUM,ASUM, EPS ,DT,DX( 50 ) ,DMAX(35 ) 

DIMENSION  ZB(100) 

C 

SSUM=0 . 

DO  100  M=1,MM 
ZB(M)=ZBL(M) 

SSUM=SSUM+ZB(M) 

100  CONTINUE 

IF(SSUM.LE.O.O)  GO  TO  105 
DO  104  M=1 ,MM 
PB1(K,M)=ZB(M)/SSUM 

104  CONTINUE 
RETURN 

105  CONTINUE 

DO  106  M=1,MM 

PB1(K,M)=PB0(K,M) 

106  CONTINUE 
RETURN 
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SUBROUTINE  TRANSP  ( K , J , BPER , SLP ,QE , AREA, ALP , BET , DC , DCOH  ) 

THIS  SUBROUTINE  IS  USED  TO  DETERMINE  THE  SEDIMENT  TRANSPORT 
CAPACTI Y  IN  THE  CHANNEL  BASED  ON  THE  MEYER-PETER,  MULLER 
BED  LOAD  EQUATION  AND  EINSTEIN’S  SUSPENDED  LOAD  EQUATION. 
TRANSPORT  CAPACITY  IS  DETERMINED  FOR  EACH  PARTICLE  SIZE. 


COMMON 

COMMON 

1 

2 

COMMON 

COMMON 

COMMON 

COMMON 

1 

COMMON 

COMMON 

COMMON 

COMMON 


/VARI /  NSEC , NWS , NPL , NDX , I MAX , MM ,QOUT( 200 ) ,QDUM( 35 ) 

/VAR2/  S£C( 35 ) , SLEN( 35 ),SL0P1(35) ,DMB( 35 , 10 ) , 

PBO(35,10),DNB(35,10),WET  K( 35 ) , POROS( 35  ) , 

SI (35) ,SW(35),SUC(35) , VFALL( 35 , 10 ) , SAREA( 35  ) 

/VAR3/  ISEC(35),IWS(35,3),IPL(35,2),IUP(35,3) 

/VAR4/  ACB(35),BEX(35),ADF(35),DELTS(35) 

/VAR5 /  V ISCO( 35 ) ,T( 35  ) 

/VAR7 /  ECB(35),GBC(35),ZBL(35),GBUP(35) ,CBLAT( 35 ) , BLAT( 35  ) , 
BTEM( 35),SB(35,5,IO),RB(35,5,10),SLOP(35,5), DDEP( 5 ) 
/VAR8/  A(35,L0),Q(I00),QP(35,10),PBl(35,10),GBO(35,10) 
/VAR9/  PORB( 35 ) ,ZSUM, ASUM , EPS , DT, DX( 50 ) ,DMAX( 35 ) 

/ VAR  1 0 /  CPR, EPR, ALAT,QLAT, Ai (35),B1(35),A2(35),B2(35) 
/VAR18/  PLAS 1(35), COHM( 35 ) 


C  DETERMINATION  OF  FLOW  CONDITIONS,  SUCH  AS  HYDRAULIC  DEPTH,  MEAN 
C  VELOCITY,  AND  BOUNDARY  SHEAR  STRESS 
AR  E  A = A  L  P*Q  E  **  BET 
WE  P  ER=C  PR* AREA** E  PR 
HYRAD=AREA/WEPER 
BPER=A2(K)*AREA**B2(K) 

VMEAN=QE/AREA 

TAU=62.4*HYRAD*SLP 

C 

I F ( PLAS I(K).GT.O.O)  THEN 

TCOH =0.0022 "PLAS I ( K ) **0 . 82 
ERATE=COHM(K)*(TAU/TCOH-l. )/ (62.4*2. 65) 

I F( ERATE.LT. 0.0)  ERATE=0 .0 
DCOH=ERATE*DT*60. 

ENDIF 

SV=SQRT(TAU/ 1.9379) 

BMV=2 . 5+VMEAN/ SV 

DC=0 . 0 
ICHECK=0 

FMAX=1 .5/(I.69+2.*ALOGI0(2.*HYRAD/DMB(K,MM)))**2 
I F ( FMAX . GT . 0 . 10)  FMAX=0.10 
IF(FMAX.LT.O.Ol)  FMAX=0.01 
FMIN=0 . 5Q*FMAX 
124  DO  204  M=1,MM 

F=l. 5/(1. 69+2 .*ALOG10(2 .*HYRAD/DMB(K,M) ) )**2 

IF(F.GT.FMAX)F=FMAX 

IF(F.LT.FMIN)F=FMIN 

TAO= . 24  2  5*F*VMEAN*VMEAN 

TTEM=TAO- 102. 96*DELTS ( K ) *DMB ( K , M ) 

C  BED  MATERIAL  LOAD  COMPUTATION 
IF(TTEM. LE.O .  )  GO  TO  127 
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1210  CONTINUE 

C  DETERMINATION  OF  RATIO  OF  SUSPENDED  BED  MATERIAL  LOAD 
FVB=VFALL(K,M) 

ZR=FVB/(0.4*SV) 

AR=2 . 0*DMAX ( K ) / HYRAD 
SUPMAX=1 . / AR 

IF(  AR  .GT.  0.9)  GO  TO  125 

CALL  POWER  (ZR,AR,FJ,SJ,1.0E-2) 
P=AR**(ZR-1.)/(11.6*(1.-AR)**ZR) 

SUSP=P*( BMV*F J+2 . 5*S J ) 

IF  (SUSP.LT.O. )  SUSP=0. 

IF(SUSP.GT.SUPMAX)  SUSP=SUPMAX 
GO  TO  126 

125  SUSP=0. 

C  DETERMINATION  OF  TRANSPORTING  CAPACITY  FOR  BED  MATERIAL  LOAD 

126  CBC(M)=( 1 .+SUSP)*WEPER*AGB(K)*TTEM**BEX(K) 
RB(K,J,M)=CBC(M)*PB1(K,M)/QE 
TTMPRE=TTEM 

M1=M 

IF(M.GE.MM)DC=1.1*DMB(K,M) 

GO  TO  204 

127  RB(K,J,M)=0. 
gbc(m)=o .0 

ICHECK=ICHECK+1 
IF(ICHECK.GT.l)  GO  TO  204 
IF(M.LE.l)  GO  TO  204 

DC=DMB ( K , M 1 ) +TTMPRE* ( DMB ( K , M ) -DMB ( K , M 1 ) ) / ( TTMPRE-TTEM ) 

204  CONTINUE 
RETURN 
END 
C 
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c 

SUBROUTINE  POWER  (Z,A,XJI ,XJ2 ,CONV) 

C 

C  THIS  SUBROUTINE  EVALUATE  J1  AND  J2  INTEGRALS 
C  USED  IN  THE  EINSTEIN  SUSPENDED  LOAD  EQUATION 
C  NOTATIONS: 

C  XJ1  =  VALUE  OF  J1  INTEGRAL 

C  XJ2  =  VALUE  OF  J2  INTEGRAL 

C  N  =  ORDER  OF  APPROXIMATION  +  1 
C  CONV  =  CONVERGENCE  CRITERION 

C 

N=1 

XJ1=0. 

XJ2=0 . 

ALG=ALOG( A) 

C=1 . 

D=-Z 
E-D+ i . 

FN=1. 

AEX=A**E 
GO  TO  102 

101  N=N+1 
C=C"'D/FN 
D=E 
E=D+1 . 

FN=FLOAT(N) 

AEX=A"*"'E 

102  IF  (ABS(E).LE. 0.001)  GO  TO  103 
XJ1=XJ1+C*(1.-AEX)/E 
XJ2=XJ2+C*( (AEX-1 . )/E**2-AEX*ALG/E) 

GO  TO  104 

103  XJ1=XJ1-C*ALG 
XJ2=XJ2-0 . 5*C*ALC**2 

104  IF  (N.EQ.l)  GO  TO  105 
CJ1=ABS( 1.-FJ1/XJ1) 

CJ2=ABS( 1.-FJ2/XJ2) 

IF  (CJ1.LE. CONV. AND. CJ2.LE. CONV)  RETURN 

105  FJ1=XJ1 
FJ2=XJ2 
GO  TO  101 
END 

C 
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C 

C 

C 

c 

c 

c 

c 


c 

c 


**************************************************************** 

SUBROUTINE  SROUT  ( K ,  J , DTX , B  PER , AREA , QE , DC , QCONC , DCOH ) 

THIS  IS  THE  SEDIMENT  ROUTING  SUBROUTINE.  IN  HERE,  THE 
LOAD  TRANSPORTED  INTO  A  CHANNEL  SECTION  IS  COMPARED  WITH 
THE  TRANSPORT  CAPACITY  AND  WHAT  IS  AVAILABLE  TO  MOVE. 

IF  AVAILABLE,  THOSE  SIZES  THAT  CAN  BE  TRANSPORTED  WILL  BE, 

WHILE  THOSE  THAT  CANNOT  WILL  NOT. 

COMMON  /VARI/  NSEG ,NWS , NPL , NDX , IMAX ,MM ,QOUT( 200 ) ,QDUM( 35 ) 

COMMON  /VAR2/  SEG( 35 ) , SLEN( 35 ) , SLOP1 ( 35 ) , DMB( 35 , 10 ) , 

1  PB0( 35 ,10) ,DNB( 35 , 10 ) ,WET  K(35 ) ,POROS(35 ) , 

2  SI (35) ,SW(35) , SUC ( 35 ) , VFALL( 35 , 10 ) , SAREA( 35 ) 

COMMON  /VAR4/  AGB( 35 ) , BEX( 35 ) , ADF( 35 ) , DELTS ( 35 ) 

COMMON  /VAR5/  VISCO( 35 ) ,T( 35 ) 

COMMON  /VAR7 /  EGB( 35 ) ,GBC( 35 ) , ZBL( 35 ) ,GBUP( 35 ) , GBLAT( 35 ) , BLAT( 35 ) , 
I  BTEM( 35 ), SB (35, 5, 10), RB (35, 5, 10), SLOP( 35 , 5 ) , DDE?( 5 ) 

COMMON  /VAR8/  A( 35 , 10 ) ,Q( 100 ) ,QP( 35 , 10 ) , PB 1 ( 35 , 10 ) ,GBO( 35 , 10 ) 
COMMON  /VAR9/  PORB( 35) ,ZSUM, ASUM, EPS , DT ,DX( 50 ) ,DMAX( 35 ) 

COMMON  /VAR10/  CPR, EPR , ALAT ,QLAT , A1 ( 35 ) , B1 ( 35 ) , A2 ( 35 ) , B2 ( 35 ) 

COMMON  /VAR18/  PLASI ( 35 ) , COHM( 35 ) 

DIMENSION  ZBNEW(IO) 

CHECK  THE  AVAILABILITY  OF  BED  MATERIAL  LOAD 

NARM=1 

PARM=1 .0 

SGBUP=0 . 

SRBE=0 . 

SBTEM=0 . 

S8LAT=0 . 

IF(ZSUM.LT.0.2)  GO  TO  50 
PARM=0 . 0 
DO  10  M=1 ,MM 
M1=MM-M+1 

PARM=PARM+PB1( 1 ,M1 ) 

IF(PB1(1,M1).LE.0.0)  GO  TO  10 
NARM=M1 
GO  TO  20 
10  CONTINUE 
20  CONTINUE 
50  DO  101  M=NARM,MM 
SGBUP=SGBUP+GBUP(M ) 

SRBE=SRBE+RB(K, J ,M) 

SBTEM=SBTEM+BTEM( M ) 

101  SBLAT=SBLAT+BLAT(M) 

SEGB=(SGBUP-SRBE*QE)*DTX-SRBE*AREA+SBTEM*A(K,J)+SBLAT 

IF(SEGB  .GE.  0.)  GO  TO  102 

EBTEM=SEGB+ZSUM*BPER 

IF(EBTEM  .GE.  0.)  GO  TO  102 

DEP=-ADF ( K )*EBTEM/ ( BPER*PARM ) 

IF(DC.GT.DMB(K,M))  GO  TO  103 
IF(DC,LE.DMB(K,l))CO  TO  102 
PC0=PB0(K,MM) 
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PCI  - PB 1 ( K , MM ) 

MMM=MM- 1 
DO  100  M=1,MMM 
M1=MM-M 

PC0=PC0+PB0(K,M1) 

PCi=PCl+PBl(K,Ml) 

IF(DC.LT.DMB(K,Ml ) )  GO  TO  100 

FAC=(DC-DMB(K,M1))/(DMB(K,M1+1)-DMB(K,M1)) 

PC0=PC0-PB0(K,M1 FAC 

PC1  =  PC1-PB  1  ( K ,  M 1 )  ••  FAC 

IF(PCl.CT.O.O)  GO  TO  97 

DEPMAX=ZSUM 

CO  TO  98 

97  DE  PMAX=2 . *DC / PC 1 

98  IF(DEPMAX.LT. ZSUM)  GO  TO  102 
IF(PCO.LE.O.O)  CO  TO  103 
DEPTEM=(2.*DC-ZSUM*PC1)/PC0 
CO  TO  99 

100  CONTINUE 
CO  TO  102 

99  IF(DEP.CT.DEPTEM)  DEP=DEPTEM 
CO  TO  103 

102  DEP^O.O 

103  CONTINUE 

I F ( PLAS I ( K ) . LE . 0 . 0 )GO  TO  104 
I F ( DEP . GT . DCOH ) DEP=DCOH 

104  ZSUM=ZSUM+DEP 
SUM=ZSUM 

I F( SUM . LE . 0 . 0  )  SUM=SEGB 
I F( SUM . LE . 0 . 0 ) SUM= .005 

C  DETERMINATION  OF  AVAILABILITY  FOR  NEXT  TIME  STEP 

DO  206  M= 1 , MM 

ZBL(M)=ZBL(M)+DEP*PB0(K,M) 

206  CONTINUE 

DO  306  M=1,MM 

ZBNEW(M)=(ZBL(M)*BPER+GBUP(M)*DTX+BTEM(M)*A(K,J)+BLAT(M)) 

1 /  (  bper+(gbc(m) / sum)*( dtx+area/ qe) ) 

305  RB(K,J,M)=ZBNEW(M)*GBC(M)/(SUM*QE) 
ECB(M)=(ZBNEW(M)-ZBL(M))*BPER 
IF(ECBCM).LE.O.O)  GO  TO  306 

IF(QE . LE . 5 .0 )  GO  TO  396 

FDEP=.5*VFALL(K,M)*BPER*DX(K)/QE 

IF(FDEP.CT.l.O)  CO  TO  306 

RBN=(  1  .-FDEP)'”(CBUP(M)  +  BLAT(M)/DTX)/QCONC 

IF(RBN . LE . RB(K , J ,M) )  CO  TO  306 

RBMAX=.99*(CBUP(M)*DTX+A(K,J)*BTEM(M)+BLAT(M))/(QE*DTX+AREA) 

IF(RBN.CT.RBMAX)RBN=RBMAX 

RB(K, J,M)=RBN 

ECB(M)=(CBUP(M)-RB(K,J,M)*QE)*DTX-AREA*RB(K,J,M) 

1  +A(K,J)*BTEM(M)+BLAT(M) 

306  CONTINUE 
RETURN 
END 
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SUBROUTINE  CEASE  ( K , J , DTX , QUP , QE , AREA ) 

THIS  SUBROUTINE  IS  USED  TO  ACCOUNT  FOR  THE  FLOW  AFTER 
THE  LATERAL  INFLOW  HAS  ENDED.  SEDIMENT  FLOWS  ARE  ZEROED. 

C 

COMMON  /VAR1/  NSEG,NWS ,NPL,NDX, IMAX,MM,QOUT(200 ) ,QDUM(35 ) 

COMMON  /VAR4/  ACB( 35 ) ,BEX( 35) ,ADF( 35 ) ,DELTS(35 ) 

COMMON  /VAR7 /  ECB(35 ) ,GBC( 35 ) ,ZBL( 35 ) ,GBUP(35 ) ,GBLAT(35 ) ,BLAT(35 ) , 
1  BTEM(35 ),SB(35,3,10),RB(35,5,10) ,SLOP(35 , 5 ) ,DDEP(5 ) 

COMMON  /VAR8/  A( 35 , 10) ,Q( 100) ,QP( 35 , 10) ,PBl(35 , 10 ) ,GBO( 35 , 10) 
COMMON  /VAR9/  PORB(35) ,ZSUM,ASUM,EPS,DT,DX(50) ,DMAX(35) 

COMMON  /VAR10/  CPR , EPR , ALAT ,QLAT , Al( 35 ) , B I ( 35 ) , A2 ( 35 ) , B2 ( 35 ) 

C 

AREA=0 . 

QE=0 . 

DO  360  M=1,MM 

EGB(M)=GBUP(M)*DTX+BTEM(M)*A(K,J)+BLAT(M) 

RB(K, J,M)=0 . 

360  CONTINUE 
RETURN 
END 
C 


81 


SUBROUTINE  OUT  ( K ,QSED ,QUP , PPM , PPPM ,  IT, ITCOM1 , DTS ) 

C 

C  THIS  SUBROUTINE  HANDLES  THE  OUTPUT  OF  THE  RESULTS. 

C  ONLY  THE  RESULTS  AT  THE  LAST  COMPUTED  UNIT  ARE 

C  WRITTEN  TO  SCREEN1 1 .OUT.  THE  OTHER  RESULTS  ARE 

C  WRITTEN  TO  FI  LEI 1  (MSEDll.OUT)  AFTER  EACH  TIME  INTERVAL. 

C 

COMMON  /VAR I  /  NSEC, NWS  ,NPL,NDX ,  IMAX  ,MM,QOUT( 200 )  ,QDUM( 35  ) 

COMMON  /VAR2 /  SEC( 35 ) , SLEN( 35 ) , SLOP1 ( 35 ) , DMB( 35 , 10 ) , 

1  PBO(35,10),DNB(35,10) , WET  K( 35 ) , POROS( 35 ) , 

2  SI(35),SW(35),SUC(35), VFALL( 35 , 1 0 ) , SAREA( 35 ) 

COMMON  /VAR6/  NUM , S IN< 35 , 5 ) , TSEDO( 35 , 10 ) , QTOUT( 35 ) ,TIME( 200 ) , 

1  TDDEP( 35 ,5), I TYPE (35) 

COMMON  /VAR17 /  TSEDEP( 35 ) , VCAP( 35 ) , VITL( 35 ) , VT( 35 ) , 

1  CBUPIN( 35 ,10),QIN(35), I  FULL  (  35 ) 

DIMENSION  QSED( 10 ) , PPPM( 10 ) , SCOUR ( 35 ) 

C 

QTOUT(K)=QTOUT(K)+(QDUM(K)+QUP)*.5*DTS/43560. 

WRITE( 1 1 , 200 , ERR=1 )K, TIME (IT) ,QUP, PPM 
WRITE( 1 1 , 2 10 , ERR=1 )(PPPM(M),M=1,MM) 

210  FORMATdX,'  PPM  BY  SIZES  :  '  , 5 ( IX , F8 . 0  ) /  16X ,  5 ( IX , F8 . 0 ) ) 

200  FORMAT ( 3X,I2,5X,F8.2,5X,F10.2,5X,F10.2) 

TSPRT  =0.0 
DO  300  M=1,MM 

TSEDO(K,M)=TSEDO(K,M)+QSED(M) 

IF(K.EQ.NSEG)  TSPRT  =  TSPRT  +  TSEDO(K,M) 

300  CONTINUE 

QDUM(K)=QUP 

IF(K.NE.NSEG. OR . IT . NE. ITCOM1 )  RETURN 
WRITE( 11 ,350,ERR=1) 

350  FORMAT ( /  / ,  '  THE  FOLLOWING  ARE  SEDIMENT  YIELDS  BY  SIZE  FRACTIONS', 

1  /,'  WITH  THE  SMALLEST  SIZE  FIRST  AND  THE  LARGEST  LAST.',/) 

DO  360  J=1 , K 

WRITE (11, 351, ERR=1 )  J 

WRITE (11, 352, ERR=1 ) (TSEDO( J ,M) ,M=1 ,MM) 

IF( ITYPE( J ) . EQ. 2 )  GO  TO  360 
SCOUR(J )=0.0 
DO  359  JJ=1,NDX 

SCOUR (J )=SCOUR( J )+TDDEP( J, JJ ) /NDX 

359  CONTINUE 

360  CONTINUE 

351  FORMAT ( / , '  THE  YIELDS  FOR  UNIT  ',12,'  IN  POUNDS.') 

352  FORMAT( 5(1X,G15.7)/5(1X,G15.7)) 

WRITE (11, 370, ERR=1 ) 

370  FORMAT ( II,'  THE  FOLLOWING  RESULTS  DISPLAY  THE  AMOUNT  OF  SEDIMENT' 

1  ,/,'  DEPOSITED  OR  REMOVED  FROM  THE  CHANNELS  AND  RESERVOIRS. 

2  , /,'  IN  THE  CASE  OF  A  CHANNEL  THE  NUMBER  REPRESENTS  THE  ' 

3  ,/,'  DEPTH  OF  DEGRADATION  OR  AGGRADATION.  FOR  A  RESEVOIR  IT 

4  ,/,'  REPRESENTS  THE  VOLUME  OF  SEDIMENT  DEPOSITED  WITH  POR-' 

5  ,/,'  OSITY  TAKEN  INTO  ACCOUNT.',/) 

DO  390  J=1,K 

IF(ITYPE(J).EQ.1)  WRITE (  11,381, ERR=1 ) J , SCOUR( J) 


IF( ITYPE( J) .EQ.2 )  WRITE (11,382, ERR=1 ) J ,TSEDEP( J ) 

390  CONTINUE 

381  FORMATC IX, '  FOR  CHANNEL  NO.  ',12,'  THE  CHANGE  IN  ELEVATION 
1G 1 5 - 7 , *  FEET.’) 

382  FORMATC  FOR  RESERVOIR  NO.  ',12,'  THE  STORED  SEDIMENT  =',C15.7 
1,*  AC-FT.') 

WR1TE( 1 i ,400 , ERR=1 ) 

400  FORMAT ( / / '  THE  FOLLOWING  ARE  THE  WATER  YIELDS  FOR  EACH  UNIT. 7/) 
WRITE( 1 1 ,401 , ERR=1 ) ( J ,QTOUT( J ) , J=1 ,NSEG ) 

401  FORMATC/'  THE  TOTAL  WATER  YIELD  FOR  UNIT  ',12,'  IS:',G15.7, 

1  'AC-FT.') 

WRITE(* , 500 )  NSEG,QTOUT(NSEG) ,TSPRT 
WRITE(6 , 500 )  NSEG,QTOUT(NSEG) ,TSPRT 
500  FORMATC 10(/),'  THE  RESULTS  FOR  THE  DOWNSTREAM  OUTLET' 

1  ,'  OF  THE  WATERSHED,  UNIT  =',I3, 

2  //'  THE  TOTAL  RUNOFF  IN  ACRE  FEET:  ',C12.5, 

3  /'  THE  TOTAL  SEDIMENT  YIELD  IN  POUNDS:  ' ,G12.5, 

4  //'  THE  SEDIMENT  YIELD  BY  SIZES:', 

5  //'  SIZE  IN  MM  POUNDS'  ) 

DO  444  M=1 ,MM 

DMBIOUT  =  DNB(NSEG ,M) 

WRITE(* ,510)  DMBIOUT, TSEDO(NSEG,M) 

WRITE(6,510)  DMBIOUT, TSEDO(NSEG,M) 

444  CONTINUE 

510  FORMATC IX, G 12 . 5 ,4X,G12 . 5 ) 

WRITE(6,520) 

WRITE(*,520) 

520  FORMATC//'  FINAL  HYDROGRAPH 7 /9X, 'TIME' ,9X, 'DISCHARGE' 

1/8X, ' (MIN. ) ' ,9X, ' (CFS) '  ) 

WRITEC*, 530 )  (TIMECi ) ,QOUT( I ) , 1=1 , ITCOM1 ) 

WRITE(f'  ,530)  (TIMECI  ),QOUT(  I ),  1  =  1,  ITCOM1) 

530  FORMATC  7X,F8.2,6X,F10.4) 

RETURN 

STOP'  **********  STOPPED  IN  OUT  **********' 

END 
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APPENDIX  B: 


FORMAT  OF 


File  Name 


FILE  1 


FILE  9 


DATA  FILES 


Table  B1 

Input  Files  for  MULTSED 


Program 


Information  Contained 


MSED1  Computational  sequence  for  subwater¬ 

shed  and  plane  units,  time  increment, 
total  time  of  hydrograph,  and 
simulation  title 

All  physical  parameters  for  subwater¬ 
shed  and  plane  units. 

MSED3  Computational  sequence  for  channel 

and  reservoirs,  time  increment,  total  time 
of  hydrograph,  and  simulation  title  (the 
time  parameters  must  agree  with  FILE  1). 
All  physical  parameters  for  channel  and 
reservoir  units. 
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Table  B2 


Format  and  Contents  of  FILE  1 
Computational  Sequence,  Plane  and  Subwatershed  Units 


Line  Variable 


Number 

Column 

5  Format 

Label 

Contents 

1 

1-10 

8A10 

TITLE 

Title  of  simulation 

2 

1-10 

F10.0 

DTIM 

Time  increment 
(min) 

11-20 

F10.0 

FTIM 

Total  length  of 
simulation  (min) 

3 

1-10 

no 

NPL 

Number  of  planes  in 
the  watershed 

11-20 

no 

NWS 

Number  of  sub¬ 
watershed  units  in 
the  watershed 

(line  4  is  repeated  for  NPL  +  NWS  units) 

4 

1-10 

no 

ISEG(I) 

Segment  number  for 
the  Ith  unit  (ISEG(I) 
must  equal  I) 

11-20 

no 

ITYPE(I) 

1  =  plane  unit, 

2  =  subwatershed  unit 

21-30 

no 

IPRINT(I) 

If  results  are  to  be 
printed,  IPRINT(I)  must 
be  positive.  If 
negative,  results  are 
not  printed. 

Table  B3 


Format  and  Contents  of  FILE  1 
Physical  Parameters  for  Subwatershed  and  Plane  Units 


Line 

Number 

Columns 

Format 

Variable 

Label 

Contents 

5 

1-80 

- 

- 

Blank 

6 

1-80 

8A10 

— 

Title  of  unit  and 
any  other 
information 

7 

1-12 

3A4 

None 

Unit  identification 

13-22 

110 

IPLANE 

1  =  plane  unit 

2  =  subwatershed  unit 

(left  plane  soil  properties) 

8 

1-10 

F10.0 

WETK(l) 

Hydraulic 

conductivity  (in. /hr) 

11-20 

F10.0 

POROS(l) 

Porosity 

(decimal  fraction) 

21-30 

F10.0 

SI  (1) 

Initial  soil 
saturation  (decimal 
fraction) 

31-30 

F10.0 

SW  (1) 

Final  soil  saturation 
(decimal  fraction) 

41-50 

F10.0 

SAVE(I) 

Average  capillary 
suction  (in.) 

51-60 

F10.0 

PLASI(I) 

Plasticity  index 
(percent) 

61-70 

F10.0 

COHM(I) 

Erosion  rate  constant 
for  cohesive  soils 
(lb/ft2-sec) 

(right  plane  soil  properties) 

9  (Same  variables  as  line  4  except  for  the  right  side  of  the  subwatershed.  If  a  plane 
unit  is  being  considered,  then  this  line  is  blank.  The  array  indices  are  all  2.) 


Table  B3  (Cont'd) 


Line  Variable 

Number  Columns  Format  Label  Contents 


(left  plane  vegetative  and  cover  characteristics) 


10 


1-10 

F10.0 

CANCOV(l) 

Percent  of  area  with 
canopy  cover 

11-20 

F10.0 

VC(1) 

Potential  storage 
volume  per  area  for 
canopy  cover  (in.) 

21-30 

F10.0 

GRNCOV(l) 

Percent  of  area 
with  ground  cover 

31-40 

F10.0 

VG(1) 

Potential  ground 
cover  storage  volume 
per  area  (in.) 

41-50 

F10.0 

PIMP(l) 

Percent  of  area  which 
has  impervious  cover 

51-60 

F10.0 

SLOOSE(l) 

Weight  of  loose  soil 
on  the  plane  (lb) 
(Usually  left  blank) 

(right  plane  vegetative  and  cover  characteristics) 


11  (Same  as  line  6  except  for  the  right  side  of  the  subwatershed.  This  line  is  blank  if  a 
plane  is  being  considered.  The  array  indices  are  all  2.) 

12  1-10  F10.0  SLOPE(l)  Overland  slope  of 

the  left  side  of  a 
subwatershed  or  a 
single  plane  unit 
(decimal  fraction). 


11-20  FI  0.0 


PLENGTH(l)  Length  of  overland 
slope  for  left 
plane  (ft). 


87 


Table  B3  (Cont’d) 


Line  Variable 

Number  Columns  Format  Label  Contents 


21-30  F10.0  DPRES(l)  Fraction  of  left 

plane  which  dees  not 
contribute  to  flow 
because  of  depres¬ 
sions  (decimal) 


(This  line  is  blank  for  a  single  plane  unit) 


13 

1-10 

F10.0 

11-20 

F10.0 

21-30 

F10.0 

14 

1-10 

F10.0 

11-20 

F10.0 

15 

1-10 

110 

11-20 

F10.0 

21-30 

F10.0 

SLOPE(2) 

Overland  slope  of 
the  right  side  of 
a  subwatershed 
unit  (decimal  fraction) 

PLENGTH(2) 

Length  of  overland 
slope  for  the  right 
plane  (ft) 

DPRES(2) 

F-action  of  right 
plane  which  does  not 
contribute  to  flow 
because  of  depression 
(decimal) 

SLOPE(3) 

Channel  slope 
(decimal  fraction) 

PLENGTH(3) 

Channel  length  (ft) 
(Note:  This  is  the 
neighboring  channel 
unit  in  the  case  of 
a  plane.) 

NRAIN 

Number  of  rainfall 
increments 

T 

Soil  temperature  (°F) 

XN 

Manning's  n  for 
channel  (blank  if 

Chezy  C  is  used) 
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Table  B3  (Cont’d) 


Line  Variable 

Number  Columns  Format  Label  Contents 


31-40 

F10.0 

A1 

A1  in  P  =  A1*A**B1 
(Note:  if  A1  is 

0.0,  the  triangular 
approximation  is 
used.) 

41-50 

F10.0 

B1 

B1  in  above 

51-60 

F10.0 

ADW 

Maximum  resistance 
parameter  for  over¬ 
land  flow 

61-70 

F10.0 

COHM(3) 

Erosion  rate  constant 
for  cohesive  soil  in 
subwatershed  channels 
(lb/ft2-sec) 

1-10 

F10.0 

RAINOLD(I) 

Rainfall  intensity 
(in./hr) 

11-20 

F10.0 

RAINT(I) 

Ending  time  of 
rainfall  intensity 
(min) 

(This  line  is  repeated  from  1  to  NRAIN  times) 

(sediment  properties) 

1-10 

F10.0 

DCOEFF 

Rainfall  splash 

detachment 

coefficient 

11-20 

F10.0 

DOF 

Overland  flow 

detachment 

coefficient 

21-30 

F10.0 

CHDOF 

Channel  detachment 

coefficient  (blank 
if  single  plane  unit) 
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Table  B3  (Cont'd) 


Line 

Number 

Columns 

Format 

Variable 

Label 

Contents 

31-40 

F10.0 

PLASI(3) 

Plasticity  index  for 
soil  in  subwatershed 
channel  (percent) 

41-50 

no 

NSED 

Number  of  sediment 
sizes  (this  number 
must  be  identical 
for  all  units) 

(sediment  grain  size  and  percentage  listed  from  smallest  to  largest) 


18 

1-10 

F10.0 

D(I) 

Sediment  diameter 
(mm) 

11-20 

F10.0 

P(I) 

Fraction  of  sediment 
equal  to  or  finer 
than  this  size 
(decimal) 

(This  line  is  repeated  from  1  to  NSED  times.) 

The  same  order  of  inputs  in  lines  5  through  18  is  followed  for  the  next 
plane  or  subwatershed  unit. 
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Table  B4 


Format  and  Contents  of  FILE  9  Computational  Sequence, 
Channel  and  Reservoir  Units 


Line 

Number 

Columns 

Format 

Variable 

Label 

Contents 

1 

1-80 

8A10 

TITLE 

Title  of  simulation 

2 

1-10 

F10.0 

DTIM 

Time  increment 
(min) 

11-20 

F10.0 

FTIM 

Total  length  of 
simulation  (min) 

(DTIM  and  FTIM  must  agree  with  FILE  1  values) 

3 

1-10 

no 

NPL 

Number  of  planes 
(must  agree  with 

FILE  1) 

11-20 

no 

NWS 

Number  of  sub- 
watersheds  (must 
agree  with  FILE  1) 

21-30 

no 

NRES 

Number  of  reservoir 
units 

31-40 

no 

NCH 

Number  of  channel 
units 

41-50 

no 

NSED 

Number  of  sediment 
sizes  (must  agree 
with  FILE  1) 

(computational  sequence) 

4 

1-5 

15 

ISEG(I) 

Channel  or  reservoir 
identification  number 
[ISEG(I)  must  equal 

1] 

6-20 

315 

IWS(I,J) 

Identification  number 
of  all  subwatershed 
units  that  flow  into 
ISEG(I)  (maximum 
of  three) 
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Table  B4  (Coat'd) 


Line 

Number 


Variable 

Columns  Format  Label  Contents 


21-30  215  IPL(I,J)  Identification  number 

of  plane  units  which 
flow  into  ISEG(I) 

[for  the  case  when 
ISEG(I)  is  a 
reservoir  no  plane 
units  are  allowed] 

31-45  315  IUP(I,J)  Identification  number 

of  upstream  input 
units  which  are 
either  channels  or 
reservoirs  (these 
numbers  correspond  to 
the  1SEG  number 
assigned  earlier  to 
the  inflowing  unit) 


Line  4  is  repeated  for  each  of  the  channel 
and  reservoir  units  (a  total  of  NCH  +  NRES) 
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Table  B5 


Format  and  Contents  of  File  9  Physical  Parameters  for 


Channel  and  Reservoir  Units 

Line 

Number  Columns 

Format 

Variable 

Label 

Contents 

5 

1-80 

- 

- 

Blank 

6 

1-80 

Channel  or  reservoir 
name  and  other 
information  (this 
line  is  not  read) 

7 

1-10 

no 

ITYPE(I) 

1  =  channel, 

2  =  reservoir  unit 

(for  the  case  of  a  channel) 

8 

1-10 

F10.0 

SLEN(I) 

Slope  length  of 
channel  (ft) 

11-20 

F10.0 

SLOP  1(1) 

Slope  of  channel 
(decimal  fraction) 

21-30 

F10.0 

WETK(I) 

Hydraulic 
conductivity  in 
channel  (in./hr) 

31-40 

F10.0 

POROS(I) 

Porosity  of 
channel  bed 
(decimal  fraction) 

41-50 

F10.0 

SI(I) 

Initial  saturation 
of  channel 
(decimal  fraction) 

51-60 

F10.0 

SW(I) 

Final  saturation  of 
channel 

(decimal  fraction) 

61-70 

F10.0 

SUC(I) 

Capillary  suction  of 
channel  bed  (in.) 

71-80 

F10.0 

PLASI(I) 

Plasticity  index  of 
channel  bed  material 
(percent) 
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Table  B5  (Cont'd) 


Line 

Number 

Columns 

Format 

Variable 

Label 

Contents 

*****  \^^en  unit  ;s  a  reservoir  ***** 

(for  the  case  of  a  reservoir,  note  the  line  number) 

8 

1-10 

F10.0 

VCAP(I) 

Storage  capacity  of 
reservoir  (acre-ft) 

11-20 

10.0 

VITL(I) 

Initial  volume  stored 
in  reservoir  at  start 
of  storm  (acre-ft) 

21-30 

F10.0 

SAREA(l) 

Surface  area  of 
reservoir  (acres) 

31-40 

F10.9 

POROS(l) 

Porosity  of  deposited 
material 

(for  the  case  of  a  reservoir) 

9 

1-10 

F10.0 

D(I,J) 

Sediment  size  (mm) 
(Must  agree  with  size 
from  FILE  1) 

11-20 

F10.0 

TRPEF(I,J) 

Trap  efficiency  for 
given  sediment  size 
(decimal  fraction) 

This  line  is  repeated  from  1  to  NSED  times. 

Note: 

For  the  case  of  a  channel,  lines  4-7  are  used,  but  for  the  case  of  a  reservoir,  only 

4  and  5.  Lines  1-3  are  used  in  both  cases  and  have  identical  meanings. 

(for  the  case  of  a  channel) 

9 

1-10 

F10.0 

XN(I) 

Manning's  n  (blank 
if  Chezy  C  is  used) 

11-20 

F10.0 

T(I) 

Temperature  of 
channel  soil  (°F) 

21-30 

F10.0 

COHM(I) 

Erosion  rate  constant 
for  cohesive  soils  in 
channel  bed 
(lb/ft2-sec) 
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Table  B5  (Cont'd) 


Line 

Number 

Columns 

Format 

Variable 

Label 

Contents 

(for  the  case  of  a  channel) 

10 

1-10 

F  10.0 

A1(I) 

A1  in  P-Al* A**B1 

11-20 

F10.0 

Bl(l) 

B1  in  above 

21-30 

F 1 0.0 

A2(I) 

A 2  in  T  =  A2* A**B2 

31-40 

FI  0.0 

B2(I) 

B2  in  above 

41-50 

F10.0 

ADF(I) 

Detachment 
coefficient  for 
channel  bed 

(for  the  case  of  a  channel) 

11 

1-10 

F10.0 

D(I,  J) 

Sediment  size  (mm). 
(Must  agree  with 
sizes  from  TAPE  1 
(listed  from  smallest 
to  largest) 

11-20 

F10.0 

PF(!,J) 

Fraction  of  sediment 

equal  to  or  smaller 
than  given  diameter 
(decimal) 


Line  11  line  is  repeated  from  1  to  NSED  times 
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Tank -Autamot  1  va  RID  Cmnand  4*09' 
ATTN:  OkOTA-O 
ATTN:  DRSTA-SP 

TSAACOM  43120 
ATTN:  DRSTS  -  B 

US  Army  Proving  Cround* 

Aberdeen  Proving  Ground  21010 
ATTN:  NGB-AAI-E 
ATTN:  DAC-M1  /Z 
Dugway  Proving  Ground  *4022 
ATTN:  STEDP-PO 
ATTN:  STEOP-MT-L-E  (2) 
Electronic  Proving  Ground  65613 
ATTN:  LO  (TAC  MAC) 

Jefferson  Proving  Ground  47252 
ATTN:  BTEJP- tH 
ATIN:  STEJP-DEH 

Yuma  Proving  Ground  *5365 
ATTN:  STEYF-PL 
ATTN:  STEYP-Ft-P 

Chemical  Syataea  Laboratory  2.01C 
ATTN:  STEAF-Ft-E 
ATTN:  SrEAP-TE-M 

ATTN:  OAOAA-CLT-E 

Diractor,  USANES  391*0 
ATTN:  NE3ER 
ATTN :  NESEN 

National  Cua^d  20310 
ATTN:  NCB-AAO-AM 

BQDA  20310 
ATTN :  DAMO-TRS 

US  Naval  Academy  21412 
ATTN:  Dlv  of  Engr  4  Weapon* 

ATTN i  Political  Set  Dapt 

Chlaf ,  Naval  Operations  20374 
ATTN:  Library 

Tyndall  AFB,  n,  33403 
ATTN:  ATE SC/ EGA 

US  Amy  Ammunition  plants 
kTTN :  rac  Engr/Envr  Ofc 
Scranton  1*505 
Mlaaiaalppl  39529 
Crana  47522 
Michigan  4*0*9 
McAl  aster  74501 
Rlvarbank  95347 
Bolaton  37660 
ATIN :  SMC BO- Eh 
Milan  3*33* 

ATTN:  SAAMI-EH 
Indiana  47111 
ATTN!  SIC  IN- EK 

Institute  of  Defense  Analyaia 
Alexandria,  VA  22311 

NASA/NSTL  Aeaearch  Library  39529 

Dept  of  Transportation  Library  20590 

tnv  Protact ion  Agency  20460 

Tranaportation  Re search  Board  20411 

Defense  Personnel  Support  Ctr  19145 

Defense  General  Supply  Ctr  232*7 


1*5 

C*/*9 


